Changeset 1367 for trunk

Timestamp:

May 23, 2009 3:47:00 AM (12 years ago)

Author:

forrest

Message:

try and improve stability - also long double interior point

Location:

trunk/Clp/src

Files:

• CbcOrClpParam.cpp (modified) (2 diffs)
• ClpCholeskyBase.cpp (modified) (71 diffs)
• ClpCholeskyBase.hpp (modified) (3 diffs)
• ClpCholeskyDense.cpp (modified) (59 diffs)
• ClpCholeskyDense.hpp (modified) (1 diff)
• ClpHelperFunctions.cpp (modified) (1 diff)
• ClpHelperFunctions.hpp (modified) (1 diff)
• ClpInterior.cpp (modified) (36 diffs)
• ClpInterior.hpp (modified) (7 diffs)
• ClpMatrixBase.cpp (modified) (1 diff)
• ClpMatrixBase.hpp (modified) (2 diffs)
• ClpMessage.cpp (modified) (1 diff)
• ClpPackedMatrix.cpp (modified) (1 diff)
• ClpPackedMatrix.hpp (modified) (1 diff)
• ClpPredictorCorrector.cpp (modified) (198 diffs)
• ClpPredictorCorrector.hpp (modified) (2 diffs)
• ClpPrimalColumnSteepest.cpp (modified) (1 diff)
• ClpPrimalColumnSteepest.hpp (modified) (1 diff)
• ClpSimplex.cpp (modified) (2 diffs)
• ClpSimplexNonlinear.cpp (modified) (4 diffs)
• ClpSimplexPrimal.cpp (modified) (9 diffs)
• ClpSolve.cpp (modified) (6 diffs)

trunk/Clp/src/CbcOrClpParam.cpp

@@ -2868 +2868 @@
      ); 
   parameters[numberParameters-1].setIntValue(1);
+#ifdef CBC_KEEP_DEPRECATED
   parameters[numberParameters++]=
     CbcOrClpParam("strengthen","Create strengthened problem",
@@ -2875 +2876 @@
      "This creates a new problem by applying the root node cuts.  All tight constraints \
 will be in resulting problem"
-     ); 
+     );
+#endif 
   parameters[numberParameters++]=
     CbcOrClpParam("strong!Branching","Number of variables to look at in strong branching",

trunk/Clp/src/ClpCholeskyBase.cpp

@@ -86 +86 @@
 #else
   // actually long double
-  workDouble_ = (double *) ClpCopyOfArray((longWork *) rhs.workDouble_,numberRows_);
+  workDouble_ = reinterpret_cast<double *> (ClpCopyOfArray(reinterpret_cast<CoinWorkDouble *> (rhs.workDouble_),numberRows_));
 #endif
   link_ = ClpCopyOfArray(rhs.link_,numberRows_);
@@ -171 +171 @@
 #else
     // actually long double
-    workDouble_ = (double *) ClpCopyOfArray((longWork *) rhs.workDouble_,numberRows_);
+    workDouble_ = reinterpret_cast<double *> (ClpCopyOfArray(reinterpret_cast<CoinWorkDouble *> (rhs.workDouble_),numberRows_));
 #endif
     link_ = ClpCopyOfArray(rhs.link_,numberRows_);
@@ -195 +195 @@
    region1 is rows+columns, region2 is rows */
 void 
-ClpCholeskyBase::solveKKT (double * region1, double * region2, const double * diagonal,
-                           double diagonalScaleFactor)
+ClpCholeskyBase::solveKKT (CoinWorkDouble * region1, CoinWorkDouble * region2, const CoinWorkDouble * diagonal,
+                           CoinWorkDouble diagonalScaleFactor)
 {
   if (!doKKT_) {
@@ -202 +202 @@
     int numberColumns = model_->numberColumns();
     int numberTotal = numberRows_+numberColumns;
-    double * region1Save = new double[numberTotal];
+    CoinWorkDouble * region1Save = new CoinWorkDouble[numberTotal];
     for (iColumn=0;iColumn<numberTotal;iColumn++) {
       region1[iColumn] *= diagonal[iColumn];
@@ -209 +209 @@
     multiplyAdd(region1+numberColumns,numberRows_,-1.0,region2,1.0);
     model_->clpMatrix()->times(1.0,region1,region2);
-    double maximumRHS = maximumAbsElement(region2,numberRows_);
-    double scale=1.0;
-    double unscale=1.0;
+    CoinWorkDouble maximumRHS = maximumAbsElement(region2,numberRows_);
+    CoinWorkDouble scale=1.0;
+    CoinWorkDouble unscale=1.0;
     if (maximumRHS>1.0e-30) {
       if (maximumRHS<=0.5) {
-        double factor=2.0;
+        CoinWorkDouble factor=2.0;
         while (maximumRHS<=0.5) {
           maximumRHS*=factor;
@@ -220 +220 @@
         } /* endwhile */
       } else if (maximumRHS>=2.0&&maximumRHS<=COIN_DBL_MAX) {
-        double factor=0.5;
+        CoinWorkDouble factor=0.5;
         while (maximumRHS>=2.0) {
           maximumRHS*=factor;
@@ -246 +246 @@
     int numberColumns = model_->numberColumns();
     int numberTotal = numberColumns + numberRowsModel;
-    double * array = new double [numberRows_];
+    CoinWorkDouble * array = new CoinWorkDouble [numberRows_];
     CoinMemcpyN(region1,numberTotal,array);
     CoinMemcpyN(region2,numberRowsModel,array+numberTotal);
@@ -253 +253 @@
     int iRow;
     for (iRow=0;iRow<numberTotal;iRow++) { 
-      if (rowsDropped_[iRow]&&fabs(array[iRow])>1.0e-8) {
+      if (rowsDropped_[iRow]&&CoinAbs(array[iRow])>1.0e-8) {
         printf("row region1 %d dropped %g\n",iRow,array[iRow]);
       }
     }
     for (;iRow<numberRows_;iRow++) {
-      if (rowsDropped_[iRow]&&fabs(array[iRow])>1.0e-8) {
+      if (rowsDropped_[iRow]&&CoinAbs(array[iRow])>1.0e-8) {
         printf("row region2 %d dropped %g\n",iRow,array[iRow]);
       }
@@ -1237 +1237 @@
 #else
     // actually long double
-    workDouble_ = (double *) (new longWork[numberRows_]);
+    workDouble_ = reinterpret_cast<double *> (new CoinWorkDouble[numberRows_]);
 #endif
     diagonal_ = new longDouble[numberRows_];
@@ -1521 +1521 @@
 /* Factorize - filling in rowsDropped and returning number dropped */
 int 
-ClpCholeskyBase::factorize(const double * diagonal, int * rowsDropped) 
+ClpCholeskyBase::factorize(const CoinWorkDouble * diagonal, int * rowsDropped) 
 {
   const CoinBigIndex * columnStart = model_->clpMatrix()->getVectorStarts();
@@ -1533 +1533 @@
   int numberColumns=model_->clpMatrix()->getNumCols();
   //perturbation
-  longDouble perturbation=model_->diagonalPerturbation()*model_->diagonalNorm();
+  CoinWorkDouble perturbation=model_->diagonalPerturbation()*model_->diagonalNorm();
   //perturbation=perturbation*perturbation*100000000.0;
   if (perturbation>1.0) {
@@ -1540 +1540 @@
       std::cout <<"large perturbation "<<perturbation<<std::endl;
 #endif
-    perturbation=sqrt(perturbation);;
+    perturbation=CoinSqrt(perturbation);
     perturbation=1.0;
   }
@@ -1548 +1548 @@
   CoinZeroN(work,numberRows_);
   int newDropped=0;
-  double largest=1.0;
-  double smallest=COIN_DBL_MAX;
+  CoinWorkDouble largest=1.0;
+  CoinWorkDouble smallest=COIN_DBL_MAX;
   int numberDense=0;
   if (!doKKT_) {
-    const double * diagonalSlack = diagonal + numberColumns;
+    const CoinWorkDouble * diagonalSlack = diagonal + numberColumns;
     if (dense_)
       numberDense=dense_->numberRows();
@@ -1577 +1577 @@
       }
     }
-    longDouble delta2 = model_->delta(); // add delta*delta to diagonal
+    CoinWorkDouble delta2 = model_->delta(); // add delta*delta to diagonal
     delta2 *= delta2;
     // largest in initial matrix
-    double largest2=1.0e-20;
+    CoinWorkDouble largest2=1.0e-20;
     for (iRow=0;iRow<numberRows_;iRow++) {
       longDouble * put = sparseFactor_+choleskyStart_[iRow];
@@ -1597 +1597 @@
             CoinBigIndex start=columnStart[iColumn];
             CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
-            longDouble multiplier = diagonal[iColumn]*elementByRow[k];
+            CoinWorkDouble multiplier = diagonal[iColumn]*elementByRow[k];
             for (CoinBigIndex j=start;j<end;j++) {
               int jRow=row[j];
               int jNewRow = permuteInverse_[jRow];
               if (jNewRow>=iRow&&!rowsDropped_[jRow]) {
-                longDouble value=element[j]*multiplier;
+                CoinWorkDouble value=element[j]*multiplier;
                 work[jNewRow] += value;
               }
@@ -1609 +1609 @@
         }
         diagonal_[iRow]=work[iRow];
-        largest2 = CoinMax(largest2,fabs(work[iRow]));
+        largest2 = CoinMax(largest2,CoinAbs(work[iRow]));
         work[iRow]=0.0;
         int j;
@@ -1615 +1615 @@
           int jRow = which[j];
           put[j]=work[jRow];
-          largest2 = CoinMax(largest2,fabs(work[jRow]));
+          largest2 = CoinMax(largest2,CoinAbs(work[jRow]));
           work[jRow]=0.0;
         }
@@ -1629 +1629 @@
     //check sizes
     largest2*=1.0e-20;
-    largest = CoinMin(largest2,1.0e-11);
+    largest = CoinMin(largest2,CHOL_SMALL_VALUE);
     int numberDroppedBefore=0;
     for (iRow=0;iRow<numberRows_;iRow++) {
@@ -1636 +1636 @@
       rowsDropped[iRow]=dropped;
       if (!dropped) {
-        longDouble diagonal = diagonal_[iRow];
+        CoinWorkDouble diagonal = diagonal_[iRow];
         if (diagonal>largest2) {
           diagonal_[iRow]=diagonal+perturbation;
@@ -1663 +1663 @@
       for ( i=0;i<numberRows_;i++) {
         assert (diagonal_[i]>=0.0);
-        diagonal_[i]=sqrt(diagonal_[i]);
+        diagonal_[i]=CoinSqrt(diagonal_[i]);
       }
       // Update dense columns (just L)
@@ -1673 +1673 @@
             a[j]=0.0;
         }
-        solveLong(a,1);
+        for (i=0;i<numberRows_;i++) {
+          int iRow = permute_[i];
+          workDouble_[i] = a[iRow];
+        }
+        for (i=0;i<numberRows_;i++) {
+          CoinWorkDouble value=workDouble_[i];
+          CoinBigIndex offset = indexStart_[i]-choleskyStart_[i];
+          CoinBigIndex j;
+          for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) {
+            int iRow = choleskyRow_[j+offset];
+            workDouble_[iRow] -= sparseFactor_[j]*value;
+          }
+        }
+        for (i=0;i<numberRows_;i++) {
+          int iRow = permute_[i];
+          a[iRow]=workDouble_[i]*diagonal_[i];
+        }
       }
       dense_->resetRowsDropped();
@@ -1682 +1698 @@
         const longDouble * a = denseColumn_+i*numberRows_;
         // do diagonal
-        longDouble value = denseDiagonal[i];
+        CoinWorkDouble value = denseDiagonal[i];
         const longDouble * b = denseColumn_+i*numberRows_;
         for (int k=0;k<numberRows_;k++) 
@@ -1688 +1704 @@
         denseDiagonal[i]=value;
         for (int j=i+1;j<numberDense;j++) {
-          longDouble value = 0.0;
+          CoinWorkDouble value = 0.0;
           const longDouble * b = denseColumn_+j*numberRows_;
           for (int k=0;k<numberRows_;k++) 
@@ -1778 +1794 @@
     }
     if (permuted) {
-      longDouble delta2 = model_->delta(); // add delta*delta to bottom
+      CoinWorkDouble delta2 = model_->delta(); // add delta*delta to bottom
       delta2 *= delta2;
       // Need to permute - ugly
@@ -1788 +1804 @@
           if (iOriginalRow<numberColumns) {
             iColumn=iOriginalRow;
-            longDouble value = diagonal[iColumn];
-            if (fabs(value)>1.0e-100) {
+            CoinWorkDouble value = diagonal[iColumn];
+            if (CoinAbs(value)>1.0e-100) {
               value = 1.0/value;
-              largest = CoinMax(largest,fabs(value));
+              largest = CoinMax(largest,CoinAbs(value));
               diagonal_[iRow] = -value;
               CoinBigIndex start=columnStart[iColumn];
@@ -1800 +1816 @@
                 if (kRow>iRow) {
                   work[kRow]=element[j];
-                  largest = CoinMax(largest,fabs(element[j]));
+                  largest = CoinMax(largest,CoinAbs(element[j]));
                 }
               }
@@ -1807 +1823 @@
             }
           } else if (iOriginalRow<numberTotal) {
-            longDouble value = diagonal[iOriginalRow];
-            if (fabs(value)>1.0e-100) {
+            CoinWorkDouble value = diagonal[iOriginalRow];
+            if (CoinAbs(value)>1.0e-100) {
               value = 1.0/value;
-              largest = CoinMax(largest,fabs(value));
+              largest = CoinMax(largest,CoinAbs(value));
             } else {
               value = 1.0e100;
@@ -1828 +1844 @@
               if (jNewRow>iRow) {
                 work[jNewRow]=elementByRow[j];
-                largest = CoinMax(largest,fabs(elementByRow[j]));
+                largest = CoinMax(largest,CoinAbs(elementByRow[j]));
               }
             }
@@ -1857 +1873 @@
             CoinBigIndex j;
             iColumn=iOriginalRow;
-            longDouble value = diagonal[iColumn];
-            if (fabs(value)>1.0e-100) {
+            CoinWorkDouble value = diagonal[iColumn];
+            if (CoinAbs(value)>1.0e-100) {
               value = 1.0/value;
               for (j=columnQuadraticStart[iColumn];
@@ -1870 +1886 @@
                 }
               }
-              largest = CoinMax(largest,fabs(value));
+              largest = CoinMax(largest,CoinAbs(value));
               diagonal_[iRow] = -value;
               CoinBigIndex start=columnStart[iColumn];
@@ -1879 +1895 @@
                 if (kRow>iRow) {
                   work[kRow]=element[j];
-                  largest = CoinMax(largest,fabs(element[j]));
+                  largest = CoinMax(largest,CoinAbs(element[j]));
                 }
               }
@@ -1886 +1902 @@
             }
           } else if (iOriginalRow<numberTotal) {
-            longDouble value = diagonal[iOriginalRow];
-            if (fabs(value)>1.0e-100) {
+            CoinWorkDouble value = diagonal[iOriginalRow];
+            if (CoinAbs(value)>1.0e-100) {
               value = 1.0/value;
-              largest = CoinMax(largest,fabs(value));
+              largest = CoinMax(largest,CoinAbs(value));
             } else {
               value = 1.0e100;
@@ -1907 +1923 @@
               if (jNewRow>iRow) {
                 work[jNewRow]=elementByRow[j];
-                largest = CoinMax(largest,fabs(elementByRow[j]));
+                largest = CoinMax(largest,CoinAbs(elementByRow[j]));
               }
             }
@@ -1931 +1947 @@
           longDouble * put = sparseFactor_+choleskyStart_[iColumn];
           int * which = choleskyRow_+indexStart_[iColumn];
-          longDouble value = diagonal[iColumn];
-          if (fabs(value)>1.0e-100) {
+          CoinWorkDouble value = diagonal[iColumn];
+          if (CoinAbs(value)>1.0e-100) {
             value = 1.0/value;
-            largest = CoinMax(largest,fabs(value));
+            largest = CoinMax(largest,CoinAbs(value));
             diagonal_[iColumn] = -value;
             CoinBigIndex start=columnStart[iColumn];
@@ -1942 +1958 @@
               //sparseFactor_[numberElements++]=element[j];
               work[row[j]+numberTotal]=element[j];
-              largest = CoinMax(largest,fabs(element[j]));
+              largest = CoinMax(largest,CoinAbs(element[j]));
             }
           } else {
@@ -1965 +1981 @@
           int * which = choleskyRow_+indexStart_[iColumn];
           int number = choleskyStart_[iColumn+1]-choleskyStart_[iColumn];
-          longDouble value = diagonal[iColumn];
+          CoinWorkDouble value = diagonal[iColumn];
           CoinBigIndex j;
-          if (fabs(value)>1.0e-100) {
+          if (CoinAbs(value)>1.0e-100) {
             value = 1.0/value;
             for (j=columnQuadraticStart[iColumn];
@@ -1978 +1994 @@
               }
             }
-            largest = CoinMax(largest,fabs(value));
+            largest = CoinMax(largest,CoinAbs(value));
             diagonal_[iColumn] = -value;
             CoinBigIndex start=columnStart[iColumn];
@@ -1984 +2000 @@
             for (j=start;j<end;j++) {
               work[row[j]+numberTotal]=element[j];
-              largest = CoinMax(largest,fabs(element[j]));
+              largest = CoinMax(largest,CoinAbs(element[j]));
             }
             for (j=0;j<number;j++) {
@@ -2005 +2021 @@
         longDouble * put = sparseFactor_+choleskyStart_[iColumn];
         int * which = choleskyRow_+indexStart_[iColumn];
-        longDouble value = diagonal[iColumn];
-        if (fabs(value)>1.0e-100) {
+        CoinWorkDouble value = diagonal[iColumn];
+        if (CoinAbs(value)>1.0e-100) {
           value = 1.0/value;
-          largest = CoinMax(largest,fabs(value));
+          largest = CoinMax(largest,CoinAbs(value));
         } else {
           value = 1.0e100;
@@ -2023 +2039 @@
       }
       // Finish diagonal
-      longDouble delta2 = model_->delta(); // add delta*delta to bottom
+      CoinWorkDouble delta2 = model_->delta(); // add delta*delta to bottom
       delta2 *= delta2;
       for (iRow=0;iRow<numberRowsModel;iRow++) {
@@ -2037 +2053 @@
     //check sizes
     largest*=1.0e-20;
-    largest = CoinMin(largest,1.0e-11);
+    largest = CoinMin(largest,CHOL_SMALL_VALUE);
     doubleParameters_[10]=CoinMax(1.0e-20,largest);
     integerParameters_[20]=0;
@@ -2056 +2072 @@
     // Should save adjustments in ..R_
     int n1=0,n2=0;
-    double * primalR = model_->primalR();
-    double * dualR = model_->dualR();
+    CoinWorkDouble * primalR = model_->primalR();
+    CoinWorkDouble * dualR = model_->dualR();
     for (iRow=0;iRow<numberTotal;iRow++) { 
       if (rowsDropped2[iRow]) {
@@ -2093 +2109 @@
 ClpCholeskyBase::factorizePart2(int * rowsDropped) 
 {
-  longDouble largest=doubleParameters_[3];
-  longDouble smallest=doubleParameters_[4];
+  CoinWorkDouble largest=doubleParameters_[3];
+  CoinWorkDouble smallest=doubleParameters_[4];
   int numberDropped=integerParameters_[20];
   // probably done before
@@ -2159 +2175 @@
       for ( jRow=lastRow;jRow<iRow;jRow++) {
         int jCount = jRow-lastRow;
-        longDouble diagonalValue = diagonal_[jRow];
+        CoinWorkDouble diagonalValue = diagonal_[jRow];
         CoinBigIndex start=choleskyStart_[jRow];
         CoinBigIndex end=choleskyStart_[jRow+1];
@@ -2165 +2181 @@
           jCount--;
           CoinBigIndex get = choleskyStart_[kRow]+jCount;
-          longDouble a_jk = sparseFactor_[get];
-          longDouble value1 = d[kRow]*a_jk;
+          CoinWorkDouble a_jk = sparseFactor_[get];
+          CoinWorkDouble value1 = d[kRow]*a_jk;
           diagonalValue -= a_jk*value1;
           for (CoinBigIndex j=start;j<end;j++)
@@ -2222 +2238 @@
     }
     // for each column L[*,kRow] that affects L[*,iRow] 
-    longDouble diagonalValue=diagonal_[iRow];
+    CoinWorkDouble diagonalValue=diagonal_[iRow];
     int nextRow = link_[iRow];
     int kRow=0;
@@ -2234 +2250 @@
       CoinBigIndex end=choleskyStart_[kRow+1];
       assert(k<end);
-      longDouble a_ik=sparseFactor_[k++];
-      longDouble value1=d[kRow]*a_ik;
+      CoinWorkDouble a_ik=sparseFactor_[k++];
+      CoinWorkDouble value1=d[kRow]*a_ik;
       // update first
       first[kRow]=k;
@@ -2259 +2275 @@
             CoinBigIndex j=first[kkRow];
             //int iiRow = choleskyRow_[j+indexStart_[kkRow]-choleskyStart_[kkRow]];
-            longDouble a = sparseFactor_[j];
-            longDouble dValue = d[kkRow]*a;
+            CoinWorkDouble a = sparseFactor_[j];
+            CoinWorkDouble dValue = d[kkRow]*a;
             diagonalValue -= a*dValue;
             work[kkRow]=dValue;
@@ -2271 +2287 @@
           for (int i=0;i<length;i++) {
             int lRow = choleskyRow_[currentIndex++];
-            longDouble t0 = work[lRow];
+            CoinWorkDouble t0 = work[lRow];
             for (int kkRow=kRow;kkRow<last;kkRow++) {
               CoinBigIndex j = first[kkRow]+i;
@@ -2340 +2356 @@
         for (CoinBigIndex j=start;j<end;j++) {
           int jRow = choleskyRow_[j+offset];
-          longDouble value = sparseFactor_[j]-work[jRow];
+          CoinWorkDouble value = sparseFactor_[j]-work[jRow];
           work[jRow]=0.0;
           sparseFactor_[j]= diagonalValue*value;
@@ -2395 +2411 @@
       CoinBigIndex offset = indexStart_[iRow]-choleskyStart_[iRow];
       if (clique_[iRow]<2) {
-        longDouble dValue=d[iRow];
+        CoinWorkDouble dValue=d[iRow];
         for (CoinBigIndex k=start;k<end;k++) {
           int kRow = choleskyRow_[k+offset]; 
           assert(kRow>=firstDense_);
-          longDouble a_ik=sparseFactor_[k];
-          longDouble value1=dValue*a_ik;
+          CoinWorkDouble a_ik=sparseFactor_[k];
+          CoinWorkDouble value1=dValue*a_ik;
           diagonal_[kRow] -= value1*a_ik;
           CoinBigIndex base = choleskyStart_[kRow]-kRow-1;
           for (CoinBigIndex j=k+1;j<end;j++) {
             int jRow=choleskyRow_[j+offset];
-            longDouble a_jk = sparseFactor_[j];
+            CoinWorkDouble a_jk = sparseFactor_[j];
             sparseFactor_[base+jRow] -= a_jk*value1;
           }
@@ -2411 +2427 @@
       } else if (clique_[iRow]<3) {
         // do as pair
-        longDouble dValue0=d[iRow];
-        longDouble dValue1=d[iRow+1];
+        CoinWorkDouble dValue0=d[iRow];
+        CoinWorkDouble dValue1=d[iRow+1];
         int offset1 = first[iRow+1]-start;
         // skip row
@@ -2419 +2435 @@
           int kRow = choleskyRow_[k+offset]; 
           assert(kRow>=firstDense_);
-          longDouble a_ik0=sparseFactor_[k];
-          longDouble value0=dValue0*a_ik0;
-          longDouble a_ik1=sparseFactor_[k+offset1];
-          longDouble value1=dValue1*a_ik1;
+          CoinWorkDouble a_ik0=sparseFactor_[k];
+          CoinWorkDouble value0=dValue0*a_ik0;
+          CoinWorkDouble a_ik1=sparseFactor_[k+offset1];
+          CoinWorkDouble value1=dValue1*a_ik1;
           diagonal_[kRow] -= value0*a_ik0 + value1*a_ik1;
           CoinBigIndex base = choleskyStart_[kRow]-kRow-1;
           for (CoinBigIndex j=k+1;j<end;j++) {
             int jRow=choleskyRow_[j+offset];
-            longDouble a_jk0 = sparseFactor_[j];
-            longDouble a_jk1 = sparseFactor_[j+offset1];
+            CoinWorkDouble a_jk0 = sparseFactor_[j];
+            CoinWorkDouble a_jk1 = sparseFactor_[j+offset1];
             sparseFactor_[base+jRow] -= a_jk0*value0+a_jk1*value1;
           }
@@ -2441 +2457 @@
         // maybe later get fancy on big cliques and do transpose copy
         // seems only worth going to 3 on Intel
-        longDouble dValue0=d[iRow];
-        longDouble dValue1=d[iRow+1];
-        longDouble dValue2=d[iRow+2];
+        CoinWorkDouble dValue0=d[iRow];
+        CoinWorkDouble dValue1=d[iRow+1];
+        CoinWorkDouble dValue2=d[iRow+2];
         // get offsets and skip rows
         int offset1 = first[++iRow]-start;
@@ -2450 +2466 @@
           int kRow = choleskyRow_[k+offset]; 
           assert(kRow>=firstDense_);
-          double diagonalValue=diagonal_[kRow];
-          longDouble a_ik0=sparseFactor_[k];
-          longDouble value0=dValue0*a_ik0;
-          longDouble a_ik1=sparseFactor_[k+offset1];
-          longDouble value1=dValue1*a_ik1;
-          longDouble a_ik2=sparseFactor_[k+offset2];
-          longDouble value2=dValue2*a_ik2;
+          CoinWorkDouble diagonalValue=diagonal_[kRow];
+          CoinWorkDouble a_ik0=sparseFactor_[k];
+          CoinWorkDouble value0=dValue0*a_ik0;
+          CoinWorkDouble a_ik1=sparseFactor_[k+offset1];
+          CoinWorkDouble value1=dValue1*a_ik1;
+          CoinWorkDouble a_ik2=sparseFactor_[k+offset2];
+          CoinWorkDouble value2=dValue2*a_ik2;
           CoinBigIndex base = choleskyStart_[kRow]-kRow-1;
           diagonal_[kRow] = diagonalValue - value0*a_ik0 - value1*a_ik1 - value2*a_ik2;
           for (CoinBigIndex j=k+1;j<end;j++) {
             int jRow=choleskyRow_[j+offset];
-            longDouble a_jk0 = sparseFactor_[j];
-            longDouble a_jk1 = sparseFactor_[j+offset1];
-            longDouble a_jk2 = sparseFactor_[j+offset2];
+            CoinWorkDouble a_jk0 = sparseFactor_[j];
+            CoinWorkDouble a_jk1 = sparseFactor_[j+offset1];
+            CoinWorkDouble a_jk2 = sparseFactor_[j+offset2];
             sparseFactor_[base+jRow] -= a_jk0*value0+a_jk1*value1+a_jk2*value2;
           }
@@ -2472 +2488 @@
         // maybe later get fancy on big cliques and do transpose copy
         // maybe only worth going to 3 on Intel (but may have hidden registers)
-        longDouble dValue0=d[iRow];
-        longDouble dValue1=d[iRow+1];
-        longDouble dValue2=d[iRow+2];
-        longDouble dValue3=d[iRow+3];
+        CoinWorkDouble dValue0=d[iRow];
+        CoinWorkDouble dValue1=d[iRow+1];
+        CoinWorkDouble dValue2=d[iRow+2];
+        CoinWorkDouble dValue3=d[iRow+3];
         // get offsets and skip rows
         int offset1 = first[++iRow]-start;
@@ -2483 +2499 @@
           int kRow = choleskyRow_[k+offset]; 
           assert(kRow>=firstDense_);
-          double diagonalValue=diagonal_[kRow];
-          longDouble a_ik0=sparseFactor_[k];
-          longDouble value0=dValue0*a_ik0;
-          longDouble a_ik1=sparseFactor_[k+offset1];
-          longDouble value1=dValue1*a_ik1;
-          longDouble a_ik2=sparseFactor_[k+offset2];
-          longDouble value2=dValue2*a_ik2;
-          longDouble a_ik3=sparseFactor_[k+offset3];
-          longDouble value3=dValue3*a_ik3;
+          CoinWorkDouble diagonalValue=diagonal_[kRow];
+          CoinWorkDouble a_ik0=sparseFactor_[k];
+          CoinWorkDouble value0=dValue0*a_ik0;
+          CoinWorkDouble a_ik1=sparseFactor_[k+offset1];
+          CoinWorkDouble value1=dValue1*a_ik1;
+          CoinWorkDouble a_ik2=sparseFactor_[k+offset2];
+          CoinWorkDouble value2=dValue2*a_ik2;
+          CoinWorkDouble a_ik3=sparseFactor_[k+offset3];
+          CoinWorkDouble value3=dValue3*a_ik3;
           CoinBigIndex base = choleskyStart_[kRow]-kRow-1;
           diagonal_[kRow] = diagonalValue - (value0*a_ik0 + value1*a_ik1 + value2*a_ik2+value3*a_ik3);
           for (CoinBigIndex j=k+1;j<end;j++) {
             int jRow=choleskyRow_[j+offset];
-            longDouble a_jk0 = sparseFactor_[j];
-            longDouble a_jk1 = sparseFactor_[j+offset1];
-            longDouble a_jk2 = sparseFactor_[j+offset2];
-            longDouble a_jk3 = sparseFactor_[j+offset3];
+            CoinWorkDouble a_jk0 = sparseFactor_[j];
+            CoinWorkDouble a_jk1 = sparseFactor_[j+offset1];
+            CoinWorkDouble a_jk2 = sparseFactor_[j+offset2];
+            CoinWorkDouble a_jk3 = sparseFactor_[j+offset3];
             sparseFactor_[base+jRow] -= a_jk0*value0+a_jk1*value1+a_jk2*value2+a_jk3*value3;
           }
@@ -2520 +2536 @@
     // do change;
     int numberDense = dense_->numberRows();
-    double * change = new double[numberDense];
+    CoinWorkDouble * change = new CoinWorkDouble[numberDense];
     for (i=0;i<numberDense;i++) {
       const longDouble * a = denseColumn_+i*numberRows_;
-      longDouble value =0.0;
+      CoinWorkDouble value =0.0;
       for (int iRow=0;iRow<numberRows_;iRow++) 
         value += a[iRow]*region[iRow];
@@ -2529 +2545 @@
     }
     // solve
+#if CLP_LONG_CHOLESKY>0
+    dense_->solveLong(change);
+#else
     dense_->solve(change);
+#endif
     for (i=0;i<numberDense;i++) {
       const longDouble * a = denseColumn_+i*numberRows_;
-      longDouble value = change[i];
+      CoinWorkDouble value = change[i];
       for (int iRow=0;iRow<numberRows_;iRow++) 
         region[iRow] -= value*a[iRow];
@@ -2549 +2569 @@
 #ifdef CLP_DEBUG
   double * regionX=NULL;
-  if (sizeof(longWork)!=sizeof(double)&&type==3) {
+  if (sizeof(CoinWorkDouble)!=sizeof(double)&&type==3) {
     regionX=ClpCopyOfArray(region,numberRows_);
   }
 #endif
-  longWork * work = reinterpret_cast<longWork *> (workDouble_);
+  CoinWorkDouble * work = reinterpret_cast<CoinWorkDouble *> (workDouble_);
   int i;
   CoinBigIndex j;
@@ -2563 +2583 @@
   case 1:
     for (i=0;i<numberRows_;i++) {
-      longDouble value=work[i];
+      CoinWorkDouble value=work[i];
       CoinBigIndex offset = indexStart_[i]-choleskyStart_[i];
       for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) {
@@ -2578 +2598 @@
     for (i=numberRows_-1;i>=0;i--) {
       CoinBigIndex offset = indexStart_[i]-choleskyStart_[i];
-      longDouble value=work[i]*diagonal_[i];
+      CoinWorkDouble value=work[i]*diagonal_[i];
       for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) {
         int iRow = choleskyRow_[j+offset];
@@ -2591 +2611 @@
     for (i=0;i<firstDense_;i++) {
       CoinBigIndex offset = indexStart_[i]-choleskyStart_[i];
-      longDouble value=work[i];
+      CoinWorkDouble value=work[i];
       for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) {
         int iRow = choleskyRow_[j+offset];
@@ -2609 +2629 @@
 #endif
       for (i=numberRows_-1;i>=firstDense_;i--) {
-        longDouble value=work[i];
+        CoinWorkDouble value=work[i];
         int iRow = permute_[i];
         region[iRow]=value;
@@ -2616 +2636 @@
     for (i=firstDense_-1;i>=0;i--) {
       CoinBigIndex offset = indexStart_[i]-choleskyStart_[i];
-      longDouble value=work[i]*diagonal_[i];
+      CoinWorkDouble value=work[i]*diagonal_[i];
       for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) {
         int iRow = choleskyRow_[j+offset];
@@ -2634 +2654 @@
     double largestO=0.0;
     for (i=0;i<numberRows_;i++) {
-      largestO = CoinMax(largestO,fabs(regionX[i]));
+      largestO = CoinMax(largestO,CoinAbs(regionX[i]));
     }
     for (i=0;i<numberRows_;i++) {
@@ -2642 +2662 @@
     for (i=0;i<firstDense_;i++) {
       CoinBigIndex offset = indexStart_[i]-choleskyStart_[i];
-      longDouble value=work[i];
+      CoinWorkDouble value=work[i];
       for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) {
         int iRow = choleskyRow_[j+offset];
@@ -2656 +2676 @@
       dense.solveLong(work+firstDense_);
       for (i=numberRows_-1;i>=firstDense_;i--) {
-        longDouble value=work[i];
+        CoinWorkDouble value=work[i];
         int iRow = permute_[i];
         regionX[iRow]=value;
@@ -2663 +2683 @@
     for (i=firstDense_-1;i>=0;i--) {
       CoinBigIndex offset = indexStart_[i]-choleskyStart_[i];
-      longDouble value=work[i]*diagonal_[i];
+      CoinWorkDouble value=work[i]*diagonal_[i];
       for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) {
         int iRow = choleskyRow_[j+offset];
@@ -2675 +2695 @@
     double largestV=0.0;
     for (i=0;i<numberRows_;i++) {
-      largest = CoinMax(largest,fabs(region[i]-regionX[i]));
-      largestV = CoinMax(largestV,fabs(region[i]));
+      largest = CoinMax(largest,CoinAbs(region[i]-regionX[i]));
+      largestV = CoinMax(largestV,CoinAbs(region[i]));
     }
     printf("largest difference %g, largest %g, largest original %g\n",
@@ -2684 +2704 @@
 #endif
 }
-/* solve - 1 just first half, 2 just second half - 3 both.
-   If 1 and 2 then diagonal has sqrt of inverse otherwise inverse
-*/
+#if CLP_LONG_CHOLESKY
+/* Uses factorization to solve. */
 void 
-ClpCholeskyBase::solveLong(longDouble * region, int type)
+ClpCholeskyBase::solve (CoinWorkDouble * region) 
 {
+  assert (!whichDense_) ;
+  CoinWorkDouble * work = reinterpret_cast<CoinWorkDouble *> (workDouble_);
   int i;
   CoinBigIndex j;
   for (i=0;i<numberRows_;i++) {
     int iRow = permute_[i];
-    workDouble_[i] = region[iRow];
-  }
-  switch (type) {
-  case 1:
-    for (i=0;i<numberRows_;i++) {
-      longDouble value=workDouble_[i];
-      CoinBigIndex offset = indexStart_[i]-choleskyStart_[i];
-      for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) {
-        int iRow = choleskyRow_[j+offset];
-        workDouble_[iRow] -= sparseFactor_[j]*value;
-      }
-    }
-    for (i=0;i<numberRows_;i++) {
-      int iRow = permute_[i];
-      region[iRow]=workDouble_[i]*diagonal_[i];
-    }
-    break;
-  case 2:
-    for (i=numberRows_-1;i>=0;i--) {
-      CoinBigIndex offset = indexStart_[i]-choleskyStart_[i];
-      longDouble value=workDouble_[i]*diagonal_[i];
-      for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) {
-        int iRow = choleskyRow_[j+offset];
-        value -= sparseFactor_[j]*workDouble_[iRow];
-      }
-      workDouble_[i]=value;
+    work[i] = region[iRow];
+  }
+  for (i=0;i<firstDense_;i++) {
+    CoinBigIndex offset = indexStart_[i]-choleskyStart_[i];
+    CoinWorkDouble value=work[i];
+    for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) {
+      int iRow = choleskyRow_[j+offset];
+      work[iRow] -= sparseFactor_[j]*value;
+    }
+  }
+  if (firstDense_<numberRows_) {
+    // do dense
+    ClpCholeskyDense dense;
+    // just borrow space
+    int nDense = numberRows_-firstDense_;
+    dense.reserveSpace(this,nDense);
+#if CLP_LONG_CHOLESKY!=1
+    dense.solveLong(work+firstDense_);
+#else
+    dense.solveLongWork(work+firstDense_);
+#endif
+    for (i=numberRows_-1;i>=firstDense_;i--) {
+      CoinWorkDouble value=work[i];
       int iRow = permute_[i];
       region[iRow]=value;
     }
-    break;
-  case 3:
-    for (i=0;i<firstDense_;i++) {
-      CoinBigIndex offset = indexStart_[i]-choleskyStart_[i];
-      longDouble value=workDouble_[i];
-      for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) {
-        int iRow = choleskyRow_[j+offset];
-        workDouble_[iRow] -= sparseFactor_[j]*value;
-      }
-    }
-    if (firstDense_<numberRows_) {
-      // do dense
-      ClpCholeskyDense dense;
-      // just borrow space
-      int nDense = numberRows_-firstDense_;
-      dense.reserveSpace(this,nDense);
-      dense.solveLong(workDouble_+firstDense_);
-      for (i=numberRows_-1;i>=firstDense_;i--) {
-        longDouble value=workDouble_[i];
-        int iRow = permute_[i];
-        region[iRow]=value;
-      }
-    }
-    for (i=firstDense_-1;i>=0;i--) {
-      CoinBigIndex offset = indexStart_[i]-choleskyStart_[i];
-      longDouble value=workDouble_[i]*diagonal_[i];
-      for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) {
-        int iRow = choleskyRow_[j+offset];
-        value -= sparseFactor_[j]*workDouble_[iRow];
-      }
-      workDouble_[i]=value;
-      int iRow = permute_[i];
-      region[iRow]=value;
-    }
-    break;
+  }
+  for (i=firstDense_-1;i>=0;i--) {
+    CoinBigIndex offset = indexStart_[i]-choleskyStart_[i];
+    CoinWorkDouble value=work[i]*diagonal_[i];
+    for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) {
+      int iRow = choleskyRow_[j+offset];
+      value -= sparseFactor_[j]*work[iRow];
+    }
+    work[i]=value;
+    int iRow = permute_[i];
+    region[iRow]=value;
   }
 }
-
+#endif

trunk/Clp/src/ClpCholeskyBase.hpp

@@ -6 +6 @@
 #include "CoinPragma.hpp"
 #include "CoinFinite.hpp"
+//#define CLP_LONG_CHOLESKY 0
+#ifndef CLP_LONG_CHOLESKY 
 #define CLP_LONG_CHOLESKY 0
+#endif
+/* valid combinations are 
+   CLP_LONG_CHOLESKY 0 and COIN_LONG_WORK 0
+   CLP_LONG_CHOLESKY 1 and COIN_LONG_WORK 1
+   CLP_LONG_CHOLESKY 2 and COIN_LONG_WORK 1
+*/
+#if COIN_LONG_WORK==0
+#if CLP_LONG_CHOLESKY>0
+#define CHOLESKY_BAD_COMBINATION
+#endif
+#else
+#if CLP_LONG_CHOLESKY==0
+#define CHOLESKY_BAD_COMBINATION
+#endif
+#endif
+#ifdef CHOLESKY_BAD_COMBINATION
+#  warning("Bad combination of CLP_LONG_CHOLESKY and COIN_BIG_DOUBLE/COIN_LONG_WORK");
+"Bad combination of CLP_LONG_CHOLESKY and COIN_LONG_WORK"
+#endif
 #if CLP_LONG_CHOLESKY>1
 typedef long double longDouble;
-typedef long double longWork;
+#define CHOL_SMALL_VALUE 1.0e-15
 #elif CLP_LONG_CHOLESKY==1
 typedef double longDouble;
-typedef long double longWork;
+#define CHOL_SMALL_VALUE 1.0e-11
 #else
 typedef double longDouble;
-typedef double longWork;
+#define CHOL_SMALL_VALUE 1.0e-11
 #endif
 class ClpInterior;
@@ -46 +67 @@
   /** Factorize - filling in rowsDropped and returning number dropped.
       If return code negative then out of memory */
-  virtual int factorize(const double * diagonal, int * rowsDropped) ;
+  virtual int factorize(const CoinWorkDouble * diagonal, int * rowsDropped) ;
   /** Uses factorization to solve. */
   virtual void solve (double * region) ;
   /** Uses factorization to solve. - given as if KKT.
    region1 is rows+columns, region2 is rows */
-  virtual void solveKKT (double * region1, double * region2, const double * diagonal,
-                         double diagonalScaleFactor);
+  virtual void solveKKT (CoinWorkDouble * region1, CoinWorkDouble * region2, const CoinWorkDouble * diagonal,
+                         CoinWorkDouble diagonalScaleFactor);
+#if CLP_LONG_CHOLESKY
+  /** Uses factorization to solve. */
+  virtual void solve (CoinWorkDouble * region) ;
+#endif
   //@}
 
@@ -163 +188 @@
   */
   void solve(double * region, int type);
-  void solveLong(longDouble * region, int type);
   /// Forms ADAT - returns nonzero if not enough memory
   int preOrder(bool lowerTriangular, bool includeDiagonal, bool doKKT);

trunk/Clp/src/ClpCholeskyDense.cpp

@@ -6 +6 @@
 #include "CoinPragma.hpp"
 #include "CoinHelperFunctions.hpp"
+#include "ClpHelperFunctions.hpp"
 
 #include "ClpInterior.hpp"
@@ -165 +166 @@
   CoinZeroN(sparseFactor_,sizeFactor_);
   //perturbation
-  longDouble perturbation=model_->diagonalPerturbation()*model_->diagonalNorm();
+  CoinWorkDouble perturbation=model_->diagonalPerturbation()*model_->diagonalNorm();
   perturbation=perturbation*perturbation;
   if (perturbation>1.0) {
@@ -172 +173 @@
       std::cout <<"large perturbation "<<perturbation<<std::endl;
 #endif
-    perturbation=sqrt(perturbation);;
+    perturbation=CoinSqrt(perturbation);;
     perturbation=1.0;
   }
   int iRow;
   int newDropped=0;
-  double largest=1.0;
-  double smallest=COIN_DBL_MAX;
-  longDouble delta2 = model_->delta(); // add delta*delta to diagonal
+  CoinWorkDouble largest=1.0;
+  CoinWorkDouble smallest=COIN_DBL_MAX;
+  CoinWorkDouble delta2 = model_->delta(); // add delta*delta to diagonal
   delta2 *= delta2;
   if (!doKKT_) {
@@ -187 +188 @@
     const double * diagonalSlack = diagonal + numberColumns;
     // largest in initial matrix
-    double largest2=1.0e-20;
+    CoinWorkDouble largest2=1.0e-20;
     for (iRow=0;iRow<numberRows_;iRow++) {
       if (!rowsDropped_[iRow]) {
         CoinBigIndex startRow=rowStart[iRow];
         CoinBigIndex endRow=rowStart[iRow]+rowLength[iRow];
-        longDouble diagonalValue = diagonalSlack[iRow]+delta2;
+        CoinWorkDouble diagonalValue = diagonalSlack[iRow]+delta2;
         for (CoinBigIndex k=startRow;k<endRow;k++) {
           int iColumn=column[k];
           CoinBigIndex start=columnStart[iColumn];
           CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
-          longDouble multiplier = diagonal[iColumn]*elementByRow[k];
+          CoinWorkDouble multiplier = diagonal[iColumn]*elementByRow[k];
           for (CoinBigIndex j=start;j<end;j++) {
             int jRow=row[j];
             if (!rowsDropped_[jRow]) {
               if (jRow>iRow) {
-                longDouble value=element[j]*multiplier;
+                CoinWorkDouble value=element[j]*multiplier;
                 work[jRow] += value;
               } else if (jRow==iRow) {
-                longDouble value=element[j]*multiplier;
+                CoinWorkDouble value=element[j]*multiplier;
                 diagonalValue += value;
               }
@@ -212 +213 @@
         } 
         for (int j=iRow+1;j<numberRows_;j++) 
-          largest2 = CoinMax(largest2,fabs(work[j]));
+          largest2 = CoinMax(largest2,CoinAbs(work[j]));
         diagonal_[iRow]=diagonalValue;
-        largest2 = CoinMax(largest2,fabs(diagonalValue));
+        largest2 = CoinMax(largest2,CoinAbs(diagonalValue));
       } else {
         // dropped
@@ -224 +225 @@
     //check sizes
     largest2*=1.0e-20;
-    largest = CoinMin(largest2,1.0e-11);
+    largest = CoinMin(largest2,CHOL_SMALL_VALUE);
     int numberDroppedBefore=0;
     for (iRow=0;iRow<numberRows_;iRow++) {
@@ -231 +232 @@
       rowsDropped[iRow]=dropped;
       if (!dropped) {
-        longDouble diagonal = diagonal_[iRow];
+        CoinWorkDouble diagonal = diagonal_[iRow];
         if (diagonal>largest2) {
           diagonal_[iRow]=diagonal+perturbation;
@@ -289 +290 @@
     if (!quadratic) {
       for (iColumn=0;iColumn<numberColumns;iColumn++) {
-        longDouble value = diagonal[iColumn];
-        if (fabs(value)>1.0e-100) {
+        CoinWorkDouble value = diagonal[iColumn];
+        if (CoinAbs(value)>1.0e-100) {
           value = 1.0/value;
-          largest = CoinMax(largest,fabs(value));
+          largest = CoinMax(largest,CoinAbs(value));
           diagonal_[iColumn] = -value;
           CoinBigIndex start=columnStart[iColumn];
@@ -300 +301 @@
             //sparseFactor_[numberElements++]=element[j];
             work[row[j]+numberTotal]=element[j];
-            largest = CoinMax(largest,fabs(element[j]));
+            largest = CoinMax(largest,CoinAbs(element[j]));
           }
         } else {
@@ -315 +316 @@
       const double * quadraticElement = quadratic->getElements();
       for (iColumn=0;iColumn<numberColumns;iColumn++) {
-        longDouble value = diagonal[iColumn];
+        CoinWorkDouble value = diagonal[iColumn];
         CoinBigIndex j;
-        if (fabs(value)>1.0e-100) {
+        if (CoinAbs(value)>1.0e-100) {
           value = 1.0/value;
           for (j=columnQuadraticStart[iColumn];
@@ -328 +329 @@
             }
           }
-          largest = CoinMax(largest,fabs(value));
+          largest = CoinMax(largest,CoinAbs(value));
           diagonal_[iColumn] = -value;
           CoinBigIndex start=columnStart[iColumn];
@@ -334 +335 @@
           for (j=start;j<end;j++) {
             work[row[j]+numberTotal]=element[j];
-            largest = CoinMax(largest,fabs(element[j]));
+            largest = CoinMax(largest,CoinAbs(element[j]));
           }
         } else {
@@ -346 +347 @@
     // slacks
     for (iColumn=numberColumns;iColumn<numberTotal;iColumn++) {
-      longDouble value = diagonal[iColumn];
-      if (fabs(value)>1.0e-100) {
+      CoinWorkDouble value = diagonal[iColumn];
+      if (CoinAbs(value)>1.0e-100) {
         value = 1.0/value;
-        largest = CoinMax(largest,fabs(value));
+        largest = CoinMax(largest,CoinAbs(value));
       } else {
         value = 1.0e100;
@@ -364 +365 @@
     //check sizes
     largest*=1.0e-20;
-    largest = CoinMin(largest,1.0e-11);
+    largest = CoinMin(largest,CHOL_SMALL_VALUE);
     doubleParameters_[10]=CoinMax(1.0e-20,largest);
     integerParameters_[20]=0;
@@ -387 +388 @@
     // Should save adjustments in ..R_
     int n1=0,n2=0;
-    double * primalR = model_->primalR();
-    double * dualR = model_->dualR();
+    CoinWorkDouble * primalR = model_->primalR();
+    CoinWorkDouble * dualR = model_->dualR();
     for (iRow=0;iRow<numberTotal;iRow++) { 
       if (rowsDropped2[iRow]) {
@@ -429 +430 @@
   CoinMemcpyN(yy,numberRows_,diagonal_);
   int numberDropped=0;
-  double largest=0.0;
-  double smallest=COIN_DBL_MAX;
+  CoinWorkDouble largest=0.0;
+  CoinWorkDouble smallest=COIN_DBL_MAX;
   double dropValue = doubleParameters_[10];
   int firstPositive=integerParameters_[34];
@@ -441 +442 @@
     int addOffsetNow = numberRows_-1;;
     longDouble * workNow=sparseFactor_-1+iColumn;
-    double diagonalValue = diagonal_[iColumn];
+    CoinWorkDouble diagonalValue = diagonal_[iColumn];
     for (iRow=0;iRow<iColumn;iRow++) {
       double aj = *workNow;
       addOffsetNow--;
       workNow += addOffsetNow;
-      double multiplier = workDouble_[iRow];
-      diagonalValue -= aj*aj*multiplier;
+      diagonalValue -= aj*aj*workDouble_[iRow];
     }
     bool dropColumn=false;
@@ -776 +776 @@
                 longDouble * diagonal, longDouble * work, int * rowsDropped)
 {
-  longDouble largest=doubleParameters_[3];
-  longDouble smallest=doubleParameters_[4];
+  CoinWorkDouble largest=doubleParameters_[3];
+  CoinWorkDouble smallest=doubleParameters_[4];
   double dropValue = doubleParameters_[10];
   int firstPositive=integerParameters_[34];
@@ -783 +783 @@
   int numberDropped=0;
   int i, j, k;
-  longDouble t00,temp1;
+  CoinWorkDouble t00,temp1;
   longDouble * aa;
   aa = a-BLOCK;
@@ -790 +790 @@
     t00 = aa[j];
     for (k = 0; k < j; ++k) {
-      longDouble multiplier = work[k];
+      CoinWorkDouble multiplier = work[k];
       t00 -= a[j + k * BLOCK] * a[j + k * BLOCK] * multiplier;
     }
     bool dropColumn=false;
-    longDouble useT00=t00;
+    CoinWorkDouble useT00=t00;
     if (j+rowOffset<firstPositive) {
       // must be negative
@@ -831 +831 @@
         t00=aa[i];
         for (k = 0; k < j; ++k) {
-          longDouble multiplier = work[k];
+          CoinWorkDouble multiplier = work[k];
           t00 -= a[i + k * BLOCK] * a[j + k * BLOCK] * multiplier;
         }
@@ -874 +874 @@
     for (j = 0; j < BLOCK; j +=2) {
       aa+=2*BLOCK;
-      longDouble temp0 = diagonal[j];
-      longDouble temp1 = diagonal[j+1];
+      CoinWorkDouble temp0 = diagonal[j];
+      CoinWorkDouble temp1 = diagonal[j+1];
       for (int i=0;i<BLOCK;i+=2) {
-        longDouble t00=aUnder[i+j*BLOCK];
-        longDouble t10=aUnder[i+ BLOCK + j*BLOCK];
-        longDouble t01=aUnder[i+1+j*BLOCK];
-        longDouble t11=aUnder[i+1+ BLOCK + j*BLOCK];
+        CoinWorkDouble t00=aUnder[i+j*BLOCK];
+        CoinWorkDouble t10=aUnder[i+ BLOCK + j*BLOCK];
+        CoinWorkDouble t01=aUnder[i+1+j*BLOCK];
+        CoinWorkDouble t11=aUnder[i+1+ BLOCK + j*BLOCK];
         for (int k = 0; k < j; ++k) {
-          longDouble multiplier=work[k];
-          longDouble au0 = aUnder[i + k * BLOCK]*multiplier;
-          longDouble au1 = aUnder[i + 1 + k * BLOCK]*multiplier;
-          longDouble at0 = aTri[j + k * BLOCK];
-          longDouble at1 = aTri[j + 1 + k * BLOCK];
+          CoinWorkDouble multiplier=work[k];
+          CoinWorkDouble au0 = aUnder[i + k * BLOCK]*multiplier;
+          CoinWorkDouble au1 = aUnder[i + 1 + k * BLOCK]*multiplier;
+          CoinWorkDouble at0 = aTri[j + k * BLOCK];
+          CoinWorkDouble at1 = aTri[j + 1 + k * BLOCK];
           t00 -= au0 * at0;
           t10 -= au0 * at1;
@@ -893 +893 @@
         }
         t00 *= temp0;
-        longDouble at1 = aTri[j + 1 + j*BLOCK]*work[j];
+        CoinWorkDouble at1 = aTri[j + 1 + j*BLOCK]*work[j];
         t10 -= t00 * at1;
         t01 *= temp0;
@@ -908 +908 @@
     for (j = 0; j < BLOCK; j ++) {
       aa+=BLOCK;
-      longDouble temp1 = diagonal[j];
+      CoinWorkDouble temp1 = diagonal[j];
       for (int i=0;i<nUnder;i++) {
-        longDouble t00=aUnder[i+j*BLOCK];
+        CoinWorkDouble t00=aUnder[i+j*BLOCK];
         for (int k = 0; k < j; ++k) {
-          longDouble multiplier=work[k];
+          CoinWorkDouble multiplier=work[k];
           t00 -= aUnder[i + k * BLOCK] * aTri[j + k * BLOCK] * multiplier;
         }
@@ -937 +937 @@
 #endif
   int i, j, k;
-  longDouble t00;
+  CoinWorkDouble t00;
   longDouble * aa;
 #ifdef BLOCKUNROLL
@@ -944 +944 @@
     aa = aTri-2*BLOCK;
     for (j = 0; j < BLOCK; j +=2) {
-      longDouble t00,t01,t10,t11;
+      CoinWorkDouble t00,t01,t10,t11;
       aa+=2*BLOCK;
       aUnder2+=2;
@@ -951 +951 @@
       t10=aa[j+1+BLOCK];
       for (k = 0; k < BLOCK; ++k) {
-        longDouble multiplier = work[k];
-        longDouble a0 = aUnder2[k * BLOCK];
-        longDouble a1 = aUnder2[1 + k * BLOCK];
-        longDouble x0 = a0 * multiplier;
-        longDouble x1 = a1 * multiplier;
+        CoinWorkDouble multiplier = work[k];
+        CoinWorkDouble a0 = aUnder2[k * BLOCK];
+        CoinWorkDouble a1 = aUnder2[1 + k * BLOCK];
+        CoinWorkDouble x0 = a0 * multiplier;
+        CoinWorkDouble x1 = a1 * multiplier;
         t00 -= a0 * x0;
         t01 -= a1 * x0;
@@ -969 +969 @@
         t11=aa[i+1+BLOCK];
         for (k = 0; k < BLOCK; ++k) {
-          longDouble multiplier = work[k];
-          longDouble a0 = aUnder2[k * BLOCK]*multiplier;
-          longDouble a1 = aUnder2[1 + k * BLOCK]*multiplier;
+          CoinWorkDouble multiplier = work[k];
+          CoinWorkDouble a0 = aUnder2[k * BLOCK]*multiplier;
+          CoinWorkDouble a1 = aUnder2[1 + k * BLOCK]*multiplier;
           t00 -= aUnder[i + k * BLOCK] * a0;
           t01 -= aUnder[i + k * BLOCK] * a1;
@@ -991 +991 @@
         t00=aa[i];
         for (k = 0; k < BLOCK; ++k) {
-          longDouble multiplier = work[k];
+          CoinWorkDouble multiplier = work[k];
           t00 -= aUnder[i + k * BLOCK] * aUnder[j + k * BLOCK]*multiplier;
         }
@@ -1034 +1034 @@
       aa+=2*BLOCK;
       for (i=0;i< BLOCK;i+=2) {
-        longDouble t00=aa[i+0*BLOCK];
-        longDouble t10=aa[i+1*BLOCK];
-        longDouble t01=aa[i+1+0*BLOCK];
-        longDouble t11=aa[i+1+1*BLOCK];
+        CoinWorkDouble t00=aa[i+0*BLOCK];
+        CoinWorkDouble t10=aa[i+1*BLOCK];
+        CoinWorkDouble t01=aa[i+1+0*BLOCK];
+        CoinWorkDouble t11=aa[i+1+1*BLOCK];
         for (k=0;k<BLOCK;k++) {
-          longDouble multiplier = work[k];
-          longDouble a00=aUnder[i+k*BLOCK]*multiplier;
-          longDouble a01=aUnder[i+1+k*BLOCK]*multiplier;
+          CoinWorkDouble multiplier = work[k];
+          CoinWorkDouble a00=aUnder[i+k*BLOCK]*multiplier;
+          CoinWorkDouble a01=aUnder[i+1+k*BLOCK]*multiplier;
           t00 -= a00 * above[j + 0 + k * BLOCK];
           t10 -= a00 * above[j + 1 + k * BLOCK];
@@ -1057 +1057 @@
       aa+=4*BLOCK;
       for (i=0;i< BLOCK;i+=4) {
-        longDouble t00=aa[i+0+0*BLOCK];
-        longDouble t10=aa[i+0+1*BLOCK];
-        longDouble t20=aa[i+0+2*BLOCK];
-        longDouble t30=aa[i+0+3*BLOCK];
-        longDouble t01=aa[i+1+0*BLOCK];
-        longDouble t11=aa[i+1+1*BLOCK];
-        longDouble t21=aa[i+1+2*BLOCK];
-        longDouble t31=aa[i+1+3*BLOCK];
-        longDouble t02=aa[i+2+0*BLOCK];
-        longDouble t12=aa[i+2+1*BLOCK];
-        longDouble t22=aa[i+2+2*BLOCK];
-        longDouble t32=aa[i+2+3*BLOCK];
-        longDouble t03=aa[i+3+0*BLOCK];
-        longDouble t13=aa[i+3+1*BLOCK];
-        longDouble t23=aa[i+3+2*BLOCK];
-        longDouble t33=aa[i+3+3*BLOCK];
+        CoinWorkDouble t00=aa[i+0+0*BLOCK];
+        CoinWorkDouble t10=aa[i+0+1*BLOCK];
+        CoinWorkDouble t20=aa[i+0+2*BLOCK];
+        CoinWorkDouble t30=aa[i+0+3*BLOCK];
+        CoinWorkDouble t01=aa[i+1+0*BLOCK];
+        CoinWorkDouble t11=aa[i+1+1*BLOCK];
+        CoinWorkDouble t21=aa[i+1+2*BLOCK];
+        CoinWorkDouble t31=aa[i+1+3*BLOCK];
+        CoinWorkDouble t02=aa[i+2+0*BLOCK];
+        CoinWorkDouble t12=aa[i+2+1*BLOCK];
+        CoinWorkDouble t22=aa[i+2+2*BLOCK];
+        CoinWorkDouble t32=aa[i+2+3*BLOCK];
+        CoinWorkDouble t03=aa[i+3+0*BLOCK];
+        CoinWorkDouble t13=aa[i+3+1*BLOCK];
+        CoinWorkDouble t23=aa[i+3+2*BLOCK];
+        CoinWorkDouble t33=aa[i+3+3*BLOCK];
         const longDouble * COIN_RESTRICT aUnderNow = aUnder+i;
         const longDouble * COIN_RESTRICT aboveNow = above+j;
         for (k=0;k<BLOCK;k++) {
-          longDouble multiplier = work[k];
-          longDouble a00=aUnderNow[0]*multiplier;
-          longDouble a01=aUnderNow[1]*multiplier;
-          longDouble a02=aUnderNow[2]*multiplier;
-          longDouble a03=aUnderNow[3]*multiplier;
+          CoinWorkDouble multiplier = work[k];
+          CoinWorkDouble a00=aUnderNow[0]*multiplier;
+          CoinWorkDouble a01=aUnderNow[1]*multiplier;
+          CoinWorkDouble a02=aUnderNow[2]*multiplier;
+          CoinWorkDouble a03=aUnderNow[3]*multiplier;
           t00 -= a00 * aboveNow[0];
           t10 -= a00 * aboveNow[1];
@@ -1125 +1125 @@
       aa+=4*BLOCK;
       for (i=0;i< n;i+=2) {
-        longDouble t00=aa[i+0*BLOCK];
-        longDouble t10=aa[i+1*BLOCK];
-        longDouble t20=aa[i+2*BLOCK];
-        longDouble t30=aa[i+3*BLOCK];
-        longDouble t01=aa[i+1+0*BLOCK];
-        longDouble t11=aa[i+1+1*BLOCK];
-        longDouble t21=aa[i+1+2*BLOCK];
-        longDouble t31=aa[i+1+3*BLOCK];
+        CoinWorkDouble t00=aa[i+0*BLOCK];
+        CoinWorkDouble t10=aa[i+1*BLOCK];
+        CoinWorkDouble t20=aa[i+2*BLOCK];
+        CoinWorkDouble t30=aa[i+3*BLOCK];
+        CoinWorkDouble t01=aa[i+1+0*BLOCK];
+        CoinWorkDouble t11=aa[i+1+1*BLOCK];
+        CoinWorkDouble t21=aa[i+1+2*BLOCK];
+        CoinWorkDouble t31=aa[i+1+3*BLOCK];
         const longDouble * COIN_RESTRICT aUnderNow = aUnder+i;
         const longDouble * COIN_RESTRICT aboveNow = above+j;
         for (k=0;k<BLOCK;k++) {
-          longDouble multiplier = work[k];
-          longDouble a00=aUnderNow[0]*multiplier;
-          longDouble a01=aUnderNow[1]*multiplier;
+          CoinWorkDouble multiplier = work[k];
+          CoinWorkDouble a00=aUnderNow[0]*multiplier;
+          CoinWorkDouble a01=aUnderNow[1]*multiplier;
           t00 -= a00 * aboveNow[0];
           t10 -= a00 * aboveNow[1];
@@ -1160 +1160 @@
       }
       if (odd) {
-        longDouble t0=aa[n+0*BLOCK];
-        longDouble t1=aa[n+1*BLOCK];
-        longDouble t2=aa[n+2*BLOCK];
-        longDouble t3=aa[n+3*BLOCK];
-        longDouble a0;
+        CoinWorkDouble t0=aa[n+0*BLOCK];
+        CoinWorkDouble t1=aa[n+1*BLOCK];
+        CoinWorkDouble t2=aa[n+2*BLOCK];
+        CoinWorkDouble t3=aa[n+3*BLOCK];
+        CoinWorkDouble a0;
         for (k=0;k<BLOCK;k++) {
           a0=aUnder[n+k*BLOCK]*work[k];
@@ -1184 +1184 @@
     aa+=BLOCK;
     for (i=0;i< nUnder;i++) {
-      longDouble t00=aa[i+0*BLOCK];
+      CoinWorkDouble t00=aa[i+0*BLOCK];
       for (k=0;k<BLOCK;k++) {
-        longDouble a00=aUnder[i+k*BLOCK]*work[k];
+        CoinWorkDouble a00=aUnder[i+k*BLOCK]*work[k];
         t00 -= a00 * above[j + k * BLOCK];
       }
@@ -1293 +1293 @@
   if (numberRows_<200) {
     for (int i=0;i<numberRows_;i++) {
-      assert(fabs(region[i]-region2[i])<1.0e-3);
+      assert(CoinAbs(region[i]-region2[i])<1.0e-3);
     }
     delete [] region2;
@@ -1304 +1304 @@
 {
   int j, k;
-  longDouble t00;
+  CoinWorkDouble t00;
   for (j = 0; j < n; j ++) {
     t00 = region[j];
@@ -1322 +1322 @@
   if (n==BLOCK) {
     for (k = 0; k < BLOCK; k+=4) {
-      longDouble t0 = region2[0];
-      longDouble t1 = region2[1];
-      longDouble t2 = region2[2];
-      longDouble t3 = region2[3];
+      CoinWorkDouble t0 = region2[0];
+      CoinWorkDouble t1 = region2[1];
+      CoinWorkDouble t2 = region2[2];
+      CoinWorkDouble t3 = region2[3];
       t0 -= region[0] * a[0 + 0 * BLOCK];
       t1 -= region[0] * a[1 + 0 * BLOCK];
@@ -1416 +1416 @@
 #endif
     for (k = 0; k < n; ++k) {
-      longDouble t00 = region2[k];
+      CoinWorkDouble t00 = region2[k];
       for (j = 0; j < BLOCK; j ++) {
         t00 -= region[j] * a[k + j * BLOCK];
@@ -1431 +1431 @@
 {
   int j, k;
-  longDouble t00;
+  CoinWorkDouble t00;
   for (j = n-1; j >=0; j --) {
     t00 = region[j];
@@ -1449 +1449 @@
   if (n==BLOCK) {
     for (j = 0; j < BLOCK; j +=4) {
-      longDouble t0 = region[0];
-      longDouble t1 = region[1];
-      longDouble t2 = region[2];
-      longDouble t3 = region[3];
+      CoinWorkDouble t0 = region[0];
+      CoinWorkDouble t1 = region[1];
+      CoinWorkDouble t2 = region[2];
+      CoinWorkDouble t3 = region[3];
       t0 -= region2[0] * a[0 + 0*BLOCK];
       t1 -= region2[0] * a[0 + 1*BLOCK];
@@ -1544 +1544 @@
 #endif
     for (j = 0; j < BLOCK; j ++) {
-      longDouble t00 = region[j];
+      CoinWorkDouble t00 = region[j];
       for (k = 0; k < n; ++k) {
         t00 -= region2[k] * a[k + j * BLOCK];
@@ -1556 +1556 @@
 /* Uses factorization to solve. */
 void 
-ClpCholeskyDense::solveLong (longDouble * region) 
+ClpCholeskyDense::solveLong (CoinWorkDouble * region) 
 {
   int numberBlocks = (numberRows_+BLOCK-1)>>BLOCKSHIFT;
@@ -1619 +1619 @@
 // Forward part of solve 1
 void 
-ClpCholeskyDense::solveF1Long(longDouble * a,int n,longDouble * region)
+ClpCholeskyDense::solveF1Long(longDouble * a,int n,CoinWorkDouble * region)
 {
   int j, k;
-  longDouble t00;
+  CoinWorkDouble t00;
   for (j = 0; j < n; j ++) {
     t00 = region[j];
@@ -1634 +1634 @@
 // Forward part of solve 2
 void 
-ClpCholeskyDense::solveF2Long(longDouble * a,int n,longDouble * region, longDouble * region2)
+ClpCholeskyDense::solveF2Long(longDouble * a,int n,CoinWorkDouble * region, CoinWorkDouble * region2)
 {
   int j, k;
@@ -1640 +1640 @@
   if (n==BLOCK) {
     for (k = 0; k < BLOCK; k+=4) {
-      longDouble t0 = region2[0];
-      longDouble t1 = region2[1];
-      longDouble t2 = region2[2];
-      longDouble t3 = region2[3];
+      CoinWorkDouble t0 = region2[0];
+      CoinWorkDouble t1 = region2[1];
+      CoinWorkDouble t2 = region2[2];
+      CoinWorkDouble t3 = region2[3];
       t0 -= region[0] * a[0 + 0 * BLOCK];
       t1 -= region[0] * a[1 + 0 * BLOCK];
@@ -1734 +1734 @@
 #endif
     for (k = 0; k < n; ++k) {
-      longDouble t00 = region2[k];
+      CoinWorkDouble t00 = region2[k];
       for (j = 0; j < BLOCK; j ++) {
         t00 -= region[j] * a[k + j * BLOCK];
@@ -1746 +1746 @@
 // Backward part of solve 1
 void 
-ClpCholeskyDense::solveB1Long(longDouble * a,int n,longDouble * region)
+ClpCholeskyDense::solveB1Long(longDouble * a,int n,CoinWorkDouble * region)
 {
   int j, k;
-  longDouble t00;
+  CoinWorkDouble t00;
   for (j = n-1; j >=0; j --) {
     t00 = region[j];
@@ -1761 +1761 @@
 // Backward part of solve 2
 void 
-ClpCholeskyDense::solveB2Long(longDouble * a,int n,longDouble * region, longDouble * region2)
+ClpCholeskyDense::solveB2Long(longDouble * a,int n,CoinWorkDouble * region, CoinWorkDouble * region2)
 {
   int j, k;
@@ -1767 +1767 @@
   if (n==BLOCK) {
     for (j = 0; j < BLOCK; j +=4) {
-      longDouble t0 = region[0];
-      longDouble t1 = region[1];
-      longDouble t2 = region[2];
-      longDouble t3 = region[3];
+      CoinWorkDouble t0 = region[0];
+      CoinWorkDouble t1 = region[1];
+      CoinWorkDouble t2 = region[2];
+      CoinWorkDouble t3 = region[3];
       t0 -= region2[0] * a[0 + 0*BLOCK];
       t1 -= region2[0] * a[0 + 1*BLOCK];
@@ -1862 +1862 @@
 #endif
     for (j = 0; j < BLOCK; j ++) {
-      longDouble t00 = region[j];
+      CoinWorkDouble t00 = region[j];
       for (k = 0; k < n; ++k) {
         t00 -= region2[k] * a[k + j * BLOCK];
@@ -1874 +1874 @@
 /* Uses factorization to solve. */
 void 
-ClpCholeskyDense::solveLongWork (longWork * region) 
+ClpCholeskyDense::solveLongWork (CoinWorkDouble * region) 
 {
   int numberBlocks = (numberRows_+BLOCK-1)>>BLOCKSHIFT;
@@ -1937 +1937 @@
 // Forward part of solve 1
 void 
-ClpCholeskyDense::solveF1LongWork(longDouble * a,int n,longWork * region)
+ClpCholeskyDense::solveF1LongWork(longDouble * a,int n,CoinWorkDouble * region)
 {
   int j, k;
-  longWork t00;
+  CoinWorkDouble t00;
   for (j = 0; j < n; j ++) {
     t00 = region[j];
@@ -1952 +1952 @@
 // Forward part of solve 2
 void 
-ClpCholeskyDense::solveF2LongWork(longDouble * a,int n,longWork * region, longWork * region2)
+ClpCholeskyDense::solveF2LongWork(longDouble * a,int n,CoinWorkDouble * region, CoinWorkDouble * region2)
 {
   int j, k;
@@ -1958 +1958 @@
   if (n==BLOCK) {
     for (k = 0; k < BLOCK; k+=4) {
-      longWork t0 = region2[0];
-      longWork t1 = region2[1];
-      longWork t2 = region2[2];
-      longWork t3 = region2[3];
+      CoinWorkDouble t0 = region2[0];
+      CoinWorkDouble t1 = region2[1];
+      CoinWorkDouble t2 = region2[2];
+      CoinWorkDouble t3 = region2[3];
       t0 -= region[0] * a[0 + 0 * BLOCK];
       t1 -= region[0] * a[1 + 0 * BLOCK];
@@ -2052 +2052 @@
 #endif
     for (k = 0; k < n; ++k) {
-      longWork t00 = region2[k];
+      CoinWorkDouble t00 = region2[k];
       for (j = 0; j < BLOCK; j ++) {
         t00 -= region[j] * a[k + j * BLOCK];
@@ -2064 +2064 @@
 // Backward part of solve 1
 void 
-ClpCholeskyDense::solveB1LongWork(longDouble * a,int n,longWork * region)
+ClpCholeskyDense::solveB1LongWork(longDouble * a,int n,CoinWorkDouble * region)
 {
   int j, k;
-  longWork t00;
+  CoinWorkDouble t00;
   for (j = n-1; j >=0; j --) {
     t00 = region[j];
@@ -2079 +2079 @@
 // Backward part of solve 2
 void 
-ClpCholeskyDense::solveB2LongWork(longDouble * a,int n,longWork * region, longWork * region2)
+ClpCholeskyDense::solveB2LongWork(longDouble * a,int n,CoinWorkDouble * region, CoinWorkDouble * region2)
 {
   int j, k;
@@ -2085 +2085 @@
   if (n==BLOCK) {
     for (j = 0; j < BLOCK; j +=4) {
-      longWork t0 = region[0];
-      longWork t1 = region[1];
-      longWork t2 = region[2];
-      longWork t3 = region[3];
+      CoinWorkDouble t0 = region[0];
+      CoinWorkDouble t1 = region[1];
+      CoinWorkDouble t2 = region[2];
+      CoinWorkDouble t3 = region[3];
       t0 -= region2[0] * a[0 + 0*BLOCK];
       t1 -= region2[0] * a[0 + 1*BLOCK];
@@ -2180 +2180 @@
 #endif
     for (j = 0; j < BLOCK; j ++) {
-      longWork t00 = region[j];
+      CoinWorkDouble t00 = region[j];
       for (k = 0; k < n; ++k) {
         t00 -= region2[k] * a[k + j * BLOCK];

trunk/Clp/src/ClpCholeskyDense.hpp

@@ -89 +89 @@
   void solveB2(longDouble * a,int n,double * region,double * region2);
   /** Uses factorization to solve. */
-  void solveLong (longDouble * region) ;
+  void solveLong (CoinWorkDouble * region) ;
   /// Forward part of solve
-  void solveF1Long(longDouble * a,int n,longDouble * region);
-  void solveF2Long(longDouble * a,int n,longDouble * region,longDouble * region2);
+  void solveF1Long(longDouble * a,int n,CoinWorkDouble * region);
+  void solveF2Long(longDouble * a,int n,CoinWorkDouble * region,CoinWorkDouble * region2);
   /// Backward part of solve
-  void solveB1Long(longDouble * a,int n,longDouble * region);
-  void solveB2Long(longDouble * a,int n,longDouble * region,longDouble * region2);
+  void solveB1Long(longDouble * a,int n,CoinWorkDouble * region);
+  void solveB2Long(longDouble * a,int n,CoinWorkDouble * region,CoinWorkDouble * region2);
   /** Uses factorization to solve. */
-  void solveLongWork (longWork * region) ;
+  void solveLongWork (CoinWorkDouble * region) ;
   /// Forward part of solve
-  void solveF1LongWork(longDouble * a,int n,longWork * region);
-  void solveF2LongWork(longDouble * a,int n,longWork * region,longWork * region2);
+  void solveF1LongWork(longDouble * a,int n,CoinWorkDouble * region);
+  void solveF2LongWork(longDouble * a,int n,CoinWorkDouble * region,CoinWorkDouble * region2);
   /// Backward part of solve
-  void solveB1LongWork(longDouble * a,int n,longWork * region);
-  void solveB2LongWork(longDouble * a,int n,longWork * region,longWork * region2);
+  void solveB1LongWork(longDouble * a,int n,CoinWorkDouble * region);
+  void solveB2LongWork(longDouble * a,int n,CoinWorkDouble * region,CoinWorkDouble * region2);
   int bNumber(const longDouble * array,int &, int&);
   /// A

trunk/Clp/src/ClpHelperFunctions.cpp

@@ -112 +112 @@
   }
 }
+#if COIN_LONG_WORK 
+// For long double versions 
+CoinWorkDouble 
+maximumAbsElement(const CoinWorkDouble * region, int size)
+{
+  int i;
+  CoinWorkDouble maxValue=0.0;
+  for (i=0;i<size;i++) 
+    maxValue = CoinMax(maxValue,CoinAbs(region[i]));
+  return maxValue;
+}
+void 
+setElements(CoinWorkDouble * region, int size, CoinWorkDouble value)
+{
+  int i;
+  for (i=0;i<size;i++) 
+    region[i]=value;
+}
+void 
+multiplyAdd(const CoinWorkDouble * region1, int size, CoinWorkDouble multiplier1,
+                 CoinWorkDouble * region2, CoinWorkDouble multiplier2)
+{
+  int i;
+  if (multiplier1==1.0) {
+    if (multiplier2==1.0) {
+      for (i=0;i<size;i++) 
+        region2[i] = region1[i] + region2[i];
+    } else if (multiplier2==-1.0) {
+      for (i=0;i<size;i++) 
+        region2[i] = region1[i] - region2[i];
+    } else if (multiplier2==0.0) {
+      for (i=0;i<size;i++) 
+        region2[i] = region1[i] ;
+    } else {
+      for (i=0;i<size;i++) 
+        region2[i] = region1[i] + multiplier2*region2[i];
+    }
+  } else if (multiplier1==-1.0) {
+    if (multiplier2==1.0) {
+      for (i=0;i<size;i++) 
+        region2[i] = -region1[i] + region2[i];
+    } else if (multiplier2==-1.0) {
+      for (i=0;i<size;i++) 
+        region2[i] = -region1[i] - region2[i];
+    } else if (multiplier2==0.0) {
+      for (i=0;i<size;i++) 
+        region2[i] = -region1[i] ;
+    } else {
+      for (i=0;i<size;i++) 
+        region2[i] = -region1[i] + multiplier2*region2[i];
+    }
+  } else if (multiplier1==0.0) {
+    if (multiplier2==1.0) {
+      // nothing to do
+    } else if (multiplier2==-1.0) {
+      for (i=0;i<size;i++) 
+        region2[i] =  -region2[i];
+    } else if (multiplier2==0.0) {
+      for (i=0;i<size;i++) 
+        region2[i] =  0.0;
+    } else {
+      for (i=0;i<size;i++) 
+        region2[i] =  multiplier2*region2[i];
+    }
+  } else {
+    if (multiplier2==1.0) {
+      for (i=0;i<size;i++) 
+        region2[i] = multiplier1*region1[i] + region2[i];
+    } else if (multiplier2==-1.0) {
+      for (i=0;i<size;i++) 
+        region2[i] = multiplier1*region1[i] - region2[i];
+    } else if (multiplier2==0.0) {
+      for (i=0;i<size;i++) 
+        region2[i] = multiplier1*region1[i] ;
+    } else {
+      for (i=0;i<size;i++) 
+        region2[i] = multiplier1*region1[i] + multiplier2*region2[i];
+    }
+  }
+}
+CoinWorkDouble 
+innerProduct(const CoinWorkDouble * region1, int size, const CoinWorkDouble * region2)
+{
+  int i;
+  CoinWorkDouble value=0.0;
+  for (i=0;i<size;i++)
+    value += region1[i]*region2[i];
+  return value;
+}
+void 
+getNorms(const CoinWorkDouble * region, int size, CoinWorkDouble & norm1, CoinWorkDouble & norm2)
+{
+  norm1 = 0.0;
+  norm2 = 0.0;
+  int i;
+  for (i=0;i<size;i++) {
+    norm2 += region[i]*region[i];
+    norm1 = CoinMax(norm1,CoinAbs(region[i]));
+  }
+}
+#endif
 #ifdef DEBUG_MEMORY
 #include <malloc.h>

trunk/Clp/src/ClpHelperFunctions.hpp

@@ -18 +18 @@
 double innerProduct(const double * region1, int size, const double * region2);
 void getNorms(const double * region, int size, double & norm1, double & norm2);
+#if COIN_LONG_WORK 
+  // For long double versions 
+CoinWorkDouble maximumAbsElement(const CoinWorkDouble * region, int size);
+void setElements(CoinWorkDouble * region, int size, CoinWorkDouble value);
+void multiplyAdd(const CoinWorkDouble * region1, int size, CoinWorkDouble multiplier1,
+                 CoinWorkDouble * region2, CoinWorkDouble multiplier2);
+CoinWorkDouble innerProduct(const CoinWorkDouble * region1, int size, const CoinWorkDouble * region2);
+void getNorms(const CoinWorkDouble * region, int size, CoinWorkDouble & norm1, CoinWorkDouble & norm2);
+inline void
+CoinMemcpyN(const double * from, const int size, CoinWorkDouble * to)
+{
+  for (int i=0;i<size;i++)
+    to[i]=from[i];
+}
+inline void
+CoinMemcpyN(const CoinWorkDouble * from, const int size, double * to)
+{
+  for (int i=0;i<size;i++)
+    to[i]=from[i];
+}
+inline CoinWorkDouble
+CoinMax(const CoinWorkDouble x1, const double x2)
+{
+    return (x1 > x2) ? x1 : x2;
+}
+inline CoinWorkDouble
+CoinMax(double x1, const CoinWorkDouble x2)
+{
+    return (x1 > x2) ? x1 : x2;
+}
+inline CoinWorkDouble
+CoinMin(const CoinWorkDouble x1, const double x2)
+{
+    return (x1 < x2) ? x1 : x2;
+}
+inline CoinWorkDouble
+CoinMin(double x1, const CoinWorkDouble x2)
+{
+    return (x1 < x2) ? x1 : x2;
+}
+inline CoinWorkDouble CoinSqrt(CoinWorkDouble x)
+{
+  return sqrtl(x);
+}
+#else
+inline double CoinSqrt(double x)
+{
+  return sqrt(x);
+}
+#endif
 
 /// Following only included if ClpPdco defined

trunk/Clp/src/ClpInterior.cpp

@@ -109 +109 @@
   algorithm_(-1)
 {
-  memset(historyInfeasibility_,0,LENGTH_HISTORY*sizeof(double));
+  memset(historyInfeasibility_,0,LENGTH_HISTORY*sizeof(CoinWorkDouble));
   solveType_=3; // say interior based life form
   cholesky_ = new ClpCholeskyDense(); // put in placeholder
@@ -198 +198 @@
     algorithm_(-1)
 {
-  memset(historyInfeasibility_,0,LENGTH_HISTORY*sizeof(double));
+  memset(historyInfeasibility_,0,LENGTH_HISTORY*sizeof(CoinWorkDouble));
   solveType_=3; // say interior based life form
   cholesky_= new ClpCholeskyDense();
@@ -351 +351 @@
   algorithm_(-1)
 {
-  memset(historyInfeasibility_,0,LENGTH_HISTORY*sizeof(double));
+  memset(historyInfeasibility_,0,LENGTH_HISTORY*sizeof(CoinWorkDouble));
   solveType_=3; // say interior based life form
   cholesky_= new ClpCholeskyDense();
@@ -537 +537 @@
   int nTotal = numberRows_+numberColumns_;
   delete [] solution_;
-  solution_ = new double[nTotal];
+  solution_ = new CoinWorkDouble[nTotal];
   CoinMemcpyN(columnActivity_,  numberColumns_,solution_);
   CoinMemcpyN(rowActivity_,     numberRows_,solution_+numberColumns_);
   delete [] cost_;
-  cost_ = new double[nTotal];
+  cost_ = new CoinWorkDouble[nTotal];
   int i;
-  double direction = optimizationDirection_*objectiveScale_;
+  CoinWorkDouble direction = optimizationDirection_*objectiveScale_;
   // direction is actually scale out not scale in
   if (direction)
@@ -550 +550 @@
   for (i=0;i<numberColumns_;i++)
     cost_[i] = direction*obj[i];
-  memset(cost_+numberColumns_,0,numberRows_*sizeof(double));
+  memset(cost_+numberColumns_,0,numberRows_*sizeof(CoinWorkDouble));
   // do scaling if needed
   if (scalingFlag_>0&&!rowScale_) {
@@ -558 +558 @@
   delete [] lower_;
   delete [] upper_;
-  lower_ = new double[nTotal];
-  upper_ = new double[nTotal];
+  lower_ = new CoinWorkDouble[nTotal];
+  upper_ = new CoinWorkDouble[nTotal];
   rowLowerWork_ = lower_+numberColumns_;
   columnLowerWork_ = lower_;
@@ -587 +587 @@
   if(rowScale_) {
     for (i=0;i<numberColumns_;i++) {
-      double multiplier = rhsScale_/columnScale_[i];
+      CoinWorkDouble multiplier = rhsScale_/columnScale_[i];
       cost_[i] *= columnScale_[i];
       if (columnLowerWork_[i]>-1.0e50)
@@ -596 +596 @@
     }
     for (i=0;i<numberRows_;i++) {
-      double multiplier = rhsScale_*rowScale_[i];
+      CoinWorkDouble multiplier = rhsScale_*rowScale_[i];
       if (rowLowerWork_[i]>-1.0e50)
         rowLowerWork_[i] *= multiplier;
@@ -612 +612 @@
   }
   assert (!errorRegion_);
-  errorRegion_ = new double [numberRows_];
+  errorRegion_ = new CoinWorkDouble [numberRows_];
   assert (!rhsFixRegion_);
-  rhsFixRegion_ = new double [numberRows_];
+  rhsFixRegion_ = new CoinWorkDouble [numberRows_];
   assert (!deltaY_);
-  deltaY_ = new double [numberRows_];
+  deltaY_ = new CoinWorkDouble [numberRows_];
   CoinZeroN(deltaY_,numberRows_);
   assert (!upperSlack_);
-  upperSlack_ = new double [nTotal];
+  upperSlack_ = new CoinWorkDouble [nTotal];
   assert (!lowerSlack_);
-  lowerSlack_ = new double [nTotal];
+  lowerSlack_ = new CoinWorkDouble [nTotal];
   assert (!diagonal_);
-  diagonal_ = new double [nTotal];
+  diagonal_ = new CoinWorkDouble [nTotal];
   assert (!deltaX_);
-  deltaX_ = new double [nTotal];
+  deltaX_ = new CoinWorkDouble [nTotal];
   CoinZeroN(deltaX_,nTotal);
   assert (!deltaZ_);
-  deltaZ_ = new double [nTotal];
+  deltaZ_ = new CoinWorkDouble [nTotal];
   CoinZeroN(deltaZ_,nTotal);
   assert (!deltaW_);
-  deltaW_ = new double [nTotal];
+  deltaW_ = new CoinWorkDouble [nTotal];
   CoinZeroN(deltaW_,nTotal);
   assert (!deltaSU_);
-  deltaSU_ = new double [nTotal];
+  deltaSU_ = new CoinWorkDouble [nTotal];
   CoinZeroN(deltaSU_,nTotal);
   assert (!deltaSL_);
-  deltaSL_ = new double [nTotal];
+  deltaSL_ = new CoinWorkDouble [nTotal];
   CoinZeroN(deltaSL_,nTotal);
   assert (!primalR_);
@@ -643 +643 @@
   // create arrays if we are doing KKT
   if (cholesky_->type()>=20) {
-    primalR_ = new double [nTotal];
+    primalR_ = new CoinWorkDouble [nTotal];
     CoinZeroN(primalR_,nTotal);
-    dualR_ = new double [numberRows_];
+    dualR_ = new CoinWorkDouble [numberRows_];
     CoinZeroN(dualR_,numberRows_);
   }
   assert (!rhsB_);
-  rhsB_ = new double [numberRows_];
+  rhsB_ = new CoinWorkDouble [numberRows_];
   CoinZeroN(rhsB_,numberRows_);
   assert (!rhsU_);
-  rhsU_ = new double [nTotal];
+  rhsU_ = new CoinWorkDouble [nTotal];
   CoinZeroN(rhsU_,nTotal);
   assert (!rhsL_);
-  rhsL_ = new double [nTotal];
+  rhsL_ = new CoinWorkDouble [nTotal];
   CoinZeroN(rhsL_,nTotal);
   assert (!rhsZ_);
-  rhsZ_ = new double [nTotal];
+  rhsZ_ = new CoinWorkDouble [nTotal];
   CoinZeroN(rhsZ_,nTotal);
   assert (!rhsW_);
-  rhsW_ = new double [nTotal];
+  rhsW_ = new CoinWorkDouble [nTotal];
   CoinZeroN(rhsW_,nTotal);
   assert (!rhsC_);
-  rhsC_ = new double [nTotal];
+  rhsC_ = new CoinWorkDouble [nTotal];
   CoinZeroN(rhsC_,nTotal);
   assert (!workArray_);
-  workArray_ = new double [nTotal];
+  workArray_ = new CoinWorkDouble [nTotal];
   CoinZeroN(workArray_,nTotal);
   assert (!zVec_);
-  zVec_ = new double [nTotal];
+  zVec_ = new CoinWorkDouble [nTotal];
   CoinZeroN(zVec_,nTotal);
   assert (!wVec_);
-  wVec_ = new double [nTotal];
+  wVec_ = new CoinWorkDouble [nTotal];
   CoinZeroN(wVec_,nTotal);
   assert (!dj_);
-  dj_ = new double [nTotal];
+  dj_ = new CoinWorkDouble [nTotal];
   if (!status_)
     status_ = new unsigned char [numberRows_+numberColumns_];
@@ -687 +687 @@
   int i;
   if (optimizationDirection_!=1.0||objectiveScale_!=1.0) {
-    double scaleC = optimizationDirection_/objectiveScale_;
+    CoinWorkDouble scaleC = optimizationDirection_/objectiveScale_;
     // and modify all dual signs
     for (i=0;i<numberColumns_;i++) 
@@ -695 +695 @@
   }
   if (rowScale_) {
-    double scaleR = 1.0/rhsScale_;
+    CoinWorkDouble scaleR = 1.0/rhsScale_;
     for (i=0;i<numberColumns_;i++) {
-      double scaleFactor = columnScale_[i];
-      double valueScaled = columnActivity_[i];
+      CoinWorkDouble scaleFactor = columnScale_[i];
+      CoinWorkDouble valueScaled = columnActivity_[i];
       columnActivity_[i] = valueScaled*scaleFactor*scaleR;
-      double valueScaledDual = reducedCost_[i];
+      CoinWorkDouble valueScaledDual = reducedCost_[i];
       reducedCost_[i] = valueScaledDual/scaleFactor;
     }
     for (i=0;i<numberRows_;i++) {
-      double scaleFactor = rowScale_[i];
-      double valueScaled = rowActivity_[i];
+      CoinWorkDouble scaleFactor = rowScale_[i];
+      CoinWorkDouble valueScaled = rowActivity_[i];
       rowActivity_[i] = (valueScaled*scaleR)/scaleFactor;
-      double valueScaledDual = dual_[i];
+      CoinWorkDouble valueScaledDual = dual_[i];
       dual_[i] = valueScaledDual*scaleFactor;
     }
   } else if (rhsScale_!=1.0) {
-    double scaleR = 1.0/rhsScale_;
+    CoinWorkDouble scaleR = 1.0/rhsScale_;
     for (i=0;i<numberColumns_;i++) {
-      double valueScaled = columnActivity_[i];
+      CoinWorkDouble valueScaled = columnActivity_[i];
       columnActivity_[i] = valueScaled*scaleR;
     }
     for (i=0;i<numberRows_;i++) {
-      double valueScaled = rowActivity_[i];
+      CoinWorkDouble valueScaled = rowActivity_[i];
       rowActivity_[i] = valueScaled*scaleR;
     }
@@ -762 +762 @@
   }
   int numberBad ;
-  double largestBound, smallestBound, minimumGap;
-  double smallestObj, largestObj;
+  CoinWorkDouble largestBound, smallestBound, minimumGap;
+  CoinWorkDouble smallestObj, largestObj;
   int firstBad;
   int modifiedBounds=0;
@@ -775 +775 @@
   largestObj=0.0;
   // If bounds are too close - fix
-  double fixTolerance = 1.1*primalTolerance();
+  CoinWorkDouble fixTolerance = 1.1*primalTolerance();
   for (i=numberColumns_;i<numberColumns_+numberRows_;i++) {
-    double value;
-    value = fabs(cost_[i]);
+    CoinWorkDouble value;
+    value = CoinAbs(cost_[i]);
     if (value>1.0e50) {
       numberBad++;
@@ -805 +805 @@
     }
     if (lower_[i]>-1.0e100&&lower_[i]) {
-      value = fabs(lower_[i]);
+      value = CoinAbs(lower_[i]);
       if (value>largestBound)
         largestBound=value;
@@ -812 +812 @@
     }
     if (upper_[i]<1.0e100&&upper_[i]) {
-      value = fabs(upper_[i]);
+      value = CoinAbs(upper_[i]);
       if (value>largestBound)
         largestBound=value;
@@ -821 +821 @@
   if (largestBound)
     handler_->message(CLP_RIMSTATISTICS3,messages_)
-      <<smallestBound
-      <<largestBound
-      <<minimumGap
+      <<static_cast<double>(smallestBound)
+      <<static_cast<double>(largestBound)
+      <<static_cast<double>(minimumGap)
       <<CoinMessageEol;
   minimumGap=1.0e100;
@@ -829 +829 @@
   largestBound=0.0;
   for (i=0;i<numberColumns_;i++) {
-    double value;
-    value = fabs(cost_[i]);
+    CoinWorkDouble value;
+    value = CoinAbs(cost_[i]);
     if (value>1.0e50) {
       numberBad++;
@@ -857 +857 @@
     }
     if (lower_[i]>-1.0e100&&lower_[i]) {
-      value = fabs(lower_[i]);
+      value = CoinAbs(lower_[i]);
       if (value>largestBound)
         largestBound=value;
@@ -864 +864 @@
     }
     if (upper_[i]<1.0e100&&upper_[i]) {
-      value = fabs(upper_[i]);
+      value = CoinAbs(upper_[i]);
       if (value>largestBound)
         largestBound=value;
@@ -885 +885 @@
       <<CoinMessageEol;
   handler_->message(CLP_RIMSTATISTICS1,messages_)
-    <<smallestObj
-    <<largestObj
+    <<static_cast<double>(smallestObj)
+    <<static_cast<double>(largestObj)
     <<CoinMessageEol;  if (largestBound)
     handler_->message(CLP_RIMSTATISTICS2,messages_)
-      <<smallestBound
-      <<largestBound
-      <<minimumGap
+      <<static_cast<double>(smallestBound)
+      <<static_cast<double>(largestBound)
+      <<static_cast<double>(minimumGap)
       <<CoinMessageEol;
   return true;
@@ -991 +991 @@
 {
   int iRow,iColumn;
-  CoinMemcpyN(cost_,numberColumns_,reducedCost_);
-  matrix_->transposeTimes(-1.0,dual_,reducedCost_);
+  CoinWorkDouble * reducedCost = reinterpret_cast<CoinWorkDouble *>(reducedCost_);
+  CoinWorkDouble * dual = reinterpret_cast<CoinWorkDouble *>(dual_);
+  CoinMemcpyN(cost_,numberColumns_,reducedCost);
+  matrix_->transposeTimes(-1.0,dual,reducedCost);
   // Now modify reduced costs for quadratic
-  double quadraticOffset=quadraticDjs(reducedCost_,
+  CoinWorkDouble quadraticOffset=quadraticDjs(reducedCost,
                                       solution_,scaleFactor_);
 
@@ -1001 +1003 @@
   sumPrimalInfeasibilities_=0.0;
   sumDualInfeasibilities_=0.0;
-  double dualTolerance =  10.0*dblParam_[ClpDualTolerance];
-  double primalTolerance =  dblParam_[ClpPrimalTolerance];
-  double primalTolerance2 =  10.0*dblParam_[ClpPrimalTolerance];
+  CoinWorkDouble dualTolerance =  10.0*dblParam_[ClpDualTolerance];
+  CoinWorkDouble primalTolerance =  dblParam_[ClpPrimalTolerance];
+  CoinWorkDouble primalTolerance2 =  10.0*dblParam_[ClpPrimalTolerance];
   worstComplementarity_=0.0;
   complementarityGap_=0.0;
@@ -1009 +1011 @@
   // Done scaled - use permanent regions for output
   // but internal for bounds
-  const double * lower = lower_+numberColumns_;
-  const double * upper = upper_+numberColumns_;
+  const CoinWorkDouble * lower = lower_+numberColumns_;
+  const CoinWorkDouble * upper = upper_+numberColumns_;
   for (iRow=0;iRow<numberRows_;iRow++) {
-    double infeasibility=0.0;
-    double distanceUp = CoinMin(upper[iRow]-
-      rowActivity_[iRow],1.0e10);
-    double distanceDown = CoinMin(rowActivity_[iRow] -
-      lower[iRow],1.0e10);
+    CoinWorkDouble infeasibility=0.0;
+    CoinWorkDouble distanceUp = CoinMin(upper[iRow]-
+                                        rowActivity_[iRow],static_cast<CoinWorkDouble>(1.0e10));
+    CoinWorkDouble distanceDown = CoinMin(rowActivity_[iRow] -
+                                          lower[iRow],static_cast<CoinWorkDouble>(1.0e10));
     if (distanceUp>primalTolerance2) {
-      double value = dual_[iRow];
+      CoinWorkDouble value = dual[iRow];
       // should not be negative
       if (value<-dualTolerance) {
@@ -1029 +1031 @@
     }
     if (distanceDown>primalTolerance2) {
-      double value = dual_[iRow];
+      CoinWorkDouble value = dual[iRow];
       // should not be positive
       if (value>dualTolerance) {
@@ -1051 +1053 @@
   upper = upper_;
   for (iColumn=0;iColumn<numberColumns_;iColumn++) {
-    double infeasibility=0.0;
+    CoinWorkDouble infeasibility=0.0;
     objectiveValue_ += cost_[iColumn]*columnActivity_[iColumn];
-    double distanceUp = CoinMin(upper[iColumn]-
-      columnActivity_[iColumn],1.0e10);
-    double distanceDown = CoinMin(columnActivity_[iColumn] -
-      lower[iColumn],1.0e10);
+    CoinWorkDouble distanceUp = CoinMin(upper[iColumn]-
+                                        columnActivity_[iColumn],static_cast<CoinWorkDouble>(1.0e10));
+    CoinWorkDouble distanceDown = CoinMin(columnActivity_[iColumn] -
+                                          lower[iColumn],static_cast<CoinWorkDouble>(1.0e10));
     if (distanceUp>primalTolerance2) {
-      double value = reducedCost_[iColumn];
+      CoinWorkDouble value = reducedCost[iColumn];
       // should not be negative
       if (value<-dualTolerance) {
@@ -1069 +1071 @@
     }
     if (distanceDown>primalTolerance2) {
-      double value = reducedCost_[iColumn];
+      CoinWorkDouble value = reducedCost[iColumn];
       // should not be positive
       if (value>dualTolerance) {
@@ -1088 +1090 @@
     }
   }
+#if COIN_LONG_WORK
+  // ok as packs down
+  CoinMemcpyN(reducedCost,numberColumns_,reducedCost_);
+#endif
   // add in offset
   objectiveValue_ += 0.5*quadraticOffset;
@@ -1142 +1148 @@
 {
   // Arrays for change in columns and rhs
-  double * columnChange = new double[numberColumns_];
-  double * rowChange = new double[numberRows_];
+  CoinWorkDouble * columnChange = new CoinWorkDouble[numberColumns_];
+  CoinWorkDouble * rowChange = new CoinWorkDouble[numberRows_];
   CoinZeroN(columnChange,numberColumns_);
   CoinZeroN(rowChange,numberRows_);
   matrix_->times(1.0,columnChange,rowChange);
   int i;
-  double tolerance = primalTolerance();
+  CoinWorkDouble tolerance = primalTolerance();
   for (i=0;i<numberColumns_;i++) {
     if (columnUpper_[i]<1.0e20||columnLower_[i]>-1.0e20) { 
@@ -1154 +1160 @@
         if (fixedOrFree(i)) {
           if (columnActivity_[i]-columnLower_[i]<columnUpper_[i]-columnActivity_[i]) {
-            double change = columnLower_[i]-columnActivity_[i]; 
-            if (fabs(change)<tolerance) {
+            CoinWorkDouble change = columnLower_[i]-columnActivity_[i]; 
+            if (CoinAbs(change)<tolerance) {
               if (reallyFix)
                 columnUpper_[i]=columnLower_[i];
@@ -1162 +1168 @@
             }
           } else {
-            double change = columnUpper_[i]-columnActivity_[i]; 
-            if (fabs(change)<tolerance) {
+            CoinWorkDouble change = columnUpper_[i]-columnActivity_[i]; 
+            if (CoinAbs(change)<tolerance) {
               if (reallyFix)
                 columnLower_[i]=columnUpper_[i];
@@ -1177 +1183 @@
   matrix_->times(1.0,columnChange,rowChange);
   // If makes mess of things then don't do
-  double newSum=0.0;
+  CoinWorkDouble newSum=0.0;
   for (i=0;i<numberRows_;i++) {
-    double value=rowActivity_[i]+rowChange[i];
+    CoinWorkDouble value=rowActivity_[i]+rowChange[i];
     if (value>rowUpper_[i]+tolerance)
       newSum += value-rowUpper_[i]-tolerance;
@@ -1198 +1204 @@
             if (fixedOrFree(i+numberColumns_)) {
               if (rowActivity_[i]-rowLower_[i]<rowUpper_[i]-rowActivity_[i]) {
-                double change = rowLower_[i]-rowActivity_[i]; 
-                if (fabs(change)<tolerance) {
+                CoinWorkDouble change = rowLower_[i]-rowActivity_[i]; 
+                if (CoinAbs(change)<tolerance) {
                   if (reallyFix)
                     rowUpper_[i]=rowLower_[i];
@@ -1205 +1211 @@
                 }
               } else {
-                double change = rowLower_[i]-rowActivity_[i]; 
-                if (fabs(change)<tolerance) {
+                CoinWorkDouble change = rowLower_[i]-rowActivity_[i]; 
+                if (CoinAbs(change)<tolerance) {
                   if (reallyFix)
                     rowLower_[i]=rowUpper_[i];
@@ -1223 +1229 @@
 /* Modifies djs to allow for quadratic.
    returns quadratic offset */
-double 
-ClpInterior::quadraticDjs(double * djRegion, const double * solution,
-                          double scaleFactor)
-{
-  double quadraticOffset=0.0;
+CoinWorkDouble 
+ClpInterior::quadraticDjs(CoinWorkDouble * djRegion, const CoinWorkDouble * solution,
+                          CoinWorkDouble scaleFactor)
+{
+  CoinWorkDouble quadraticOffset=0.0;
 #ifndef NO_RTTI
   ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(objective_));
@@ -1243 +1249 @@
     int numberColumns = quadratic->getNumCols();
     for (int iColumn=0;iColumn<numberColumns;iColumn++) {
-      double value = 0.0;
+      CoinWorkDouble value = 0.0;
       for (CoinBigIndex j=columnQuadraticStart[iColumn];
            j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
         int jColumn = columnQuadratic[j];
-        double valueJ = solution[jColumn];
-        double elementValue = quadraticElement[j];
+        CoinWorkDouble valueJ = solution[jColumn];
+        CoinWorkDouble elementValue = quadraticElement[j];
         //value += valueI*valueJ*elementValue;
         value += valueJ*elementValue;

trunk/Clp/src/ClpInterior.hpp

@@ -174 +174 @@
   {algorithm_=value; } 
   /// Sum of dual infeasibilities
-  inline double sumDualInfeasibilities() const 
+  inline CoinWorkDouble sumDualInfeasibilities() const 
           { return sumDualInfeasibilities_;} 
   /// Sum of primal infeasibilities
-  inline double sumPrimalInfeasibilities() const 
+  inline CoinWorkDouble sumPrimalInfeasibilities() const 
           { return sumPrimalInfeasibilities_;} 
   /// dualObjective.
-  inline double dualObjective() const
+  inline CoinWorkDouble dualObjective() const
   { return dualObjective_;}
   /// primalObjective.
-  inline double primalObjective() const
+  inline CoinWorkDouble primalObjective() const
   { return primalObjective_;}
   /// diagonalNorm
-  inline double diagonalNorm() const
+  inline CoinWorkDouble diagonalNorm() const
   { return diagonalNorm_;}
   /// linearPerturbation
-  inline double linearPerturbation() const
+  inline CoinWorkDouble linearPerturbation() const
   { return linearPerturbation_;}
-  inline void setLinearPerturbation(double value)
+  inline void setLinearPerturbation(CoinWorkDouble value)
   { linearPerturbation_=value;}
+  /// projectionTolerance
+  inline CoinWorkDouble projectionTolerance() const
+  { return projectionTolerance_;}
+  inline void setProjectionTolerance(CoinWorkDouble value)
+  { projectionTolerance_=value;}
   /// diagonalPerturbation
-  inline double diagonalPerturbation() const
+  inline CoinWorkDouble diagonalPerturbation() const
   { return diagonalPerturbation_;}
-  inline void setDiagonalPerturbation(double value)
+  inline void setDiagonalPerturbation(CoinWorkDouble value)
   { diagonalPerturbation_=value;}
   /// gamma
-  inline double gamma() const
+  inline CoinWorkDouble gamma() const
   { return gamma_;}
-  inline void setGamma(double value)
+  inline void setGamma(CoinWorkDouble value)
   { gamma_=value;}
   /// delta
-  inline double delta() const
+  inline CoinWorkDouble delta() const
   { return delta_;}
-  inline void setDelta(double value)
+  inline void setDelta(CoinWorkDouble value)
   { delta_=value;}
   /// ComplementarityGap
-  inline double complementarityGap() const 
+  inline CoinWorkDouble complementarityGap() const 
           { return complementarityGap_;} 
   //@}
@@ -216 +221 @@
   //@{ 
   /// Largest error on Ax-b
-  inline double largestPrimalError() const
+  inline CoinWorkDouble largestPrimalError() const
           { return largestPrimalError_;} 
   /// Largest error on basic duals
-  inline double largestDualError() const
+  inline CoinWorkDouble largestDualError() const
           { return largestDualError_;} 
   /// Maximum iterations
@@ -234 +239 @@
   void fixFixed(bool reallyFix=true);
   /// Primal erturbation vector
-  inline double * primalR() const
+  inline CoinWorkDouble * primalR() const
   { return primalR_;}
   /// Dual erturbation vector
-  inline double * dualR() const
+  inline CoinWorkDouble * dualR() const
   { return dualR_;}
   //@}
@@ -260 +265 @@
   //@{
   /// Raw objective value (so always minimize)
-  inline double rawObjectiveValue() const
+  inline CoinWorkDouble rawObjectiveValue() const
   { return objectiveValue_;}
   /// Returns 1 if sequence indicates column
@@ -272 +277 @@
   /** Modifies djs to allow for quadratic.
       returns quadratic offset */
-  double quadraticDjs(double * djRegion, const double * solution,
-                      double scaleFactor);
+  CoinWorkDouble quadraticDjs(CoinWorkDouble * djRegion, const CoinWorkDouble * solution,
+                      CoinWorkDouble scaleFactor);
 
   /// To say a variable is fixed
@@ -370 +375 @@
   //@{
   /// Largest error on Ax-b
-  double largestPrimalError_;
+  CoinWorkDouble largestPrimalError_;
   /// Largest error on basic duals
-  double largestDualError_;
+  CoinWorkDouble largestDualError_;
   /// Sum of dual infeasibilities
-  double sumDualInfeasibilities_;
+  CoinWorkDouble sumDualInfeasibilities_;
   /// Sum of primal infeasibilities
-  double sumPrimalInfeasibilities_;
+  CoinWorkDouble sumPrimalInfeasibilities_;
   /// Worst complementarity
-  double worstComplementarity_;
+  CoinWorkDouble worstComplementarity_;
   /// 
 public:
-  double xsize_;
-  double zsize_;
+  CoinWorkDouble xsize_;
+  CoinWorkDouble zsize_;
 protected:
   /// Working copy of lower bounds (Owner of arrays below)
-  double * lower_;
+  CoinWorkDouble * lower_;
   /// Row lower bounds - working copy
-  double * rowLowerWork_;
+  CoinWorkDouble * rowLowerWork_;
   /// Column lower bounds - working copy
-  double * columnLowerWork_;
+  CoinWorkDouble * columnLowerWork_;
   /// Working copy of upper bounds (Owner of arrays below)
-  double * upper_;
+  CoinWorkDouble * upper_;
   /// Row upper bounds - working copy
-  double * rowUpperWork_;
+  CoinWorkDouble * rowUpperWork_;
   /// Column upper bounds - working copy
-  double * columnUpperWork_;
+  CoinWorkDouble * columnUpperWork_;
   /// Working copy of objective 
-  double * cost_;
+  CoinWorkDouble * cost_;
 public:
   /// Rhs
-  double * rhs_;
-  double * x_;
-  double * y_;
-  double * dj_;
+  CoinWorkDouble * rhs_;
+  CoinWorkDouble * x_;
+  CoinWorkDouble * y_;
+  CoinWorkDouble * dj_;
 protected:
   /// Pointer to Lsqr object
@@ -411 +416 @@
   /// Below here is standard barrier stuff
   /// mu.
-  double mu_;
+  CoinWorkDouble mu_;
   /// objectiveNorm.
-  double objectiveNorm_;
+  CoinWorkDouble objectiveNorm_;
   /// rhsNorm.
-  double rhsNorm_;
+  CoinWorkDouble rhsNorm_;
   /// solutionNorm.
-  double solutionNorm_;
+  CoinWorkDouble solutionNorm_;
   /// dualObjective.
-  double dualObjective_;
+  CoinWorkDouble dualObjective_;
   /// primalObjective.
-  double primalObjective_;
+  CoinWorkDouble primalObjective_;
   /// diagonalNorm.
-  double diagonalNorm_;
+  CoinWorkDouble diagonalNorm_;
   /// stepLength
-  double stepLength_;
+  CoinWorkDouble stepLength_;
   /// linearPerturbation
-  double linearPerturbation_;
+  CoinWorkDouble linearPerturbation_;
   /// diagonalPerturbation
-  double diagonalPerturbation_;
+  CoinWorkDouble diagonalPerturbation_;
   // gamma from Saunders and Tomlin regularized
-  double gamma_;
+  CoinWorkDouble gamma_;
   // delta from Saunders and Tomlin regularized
-  double delta_;
+  CoinWorkDouble delta_;
   /// targetGap
-  double targetGap_;
+  CoinWorkDouble targetGap_;
   /// projectionTolerance
-  double projectionTolerance_;
+  CoinWorkDouble projectionTolerance_;
   /// maximumRHSError.  maximum Ax
-  double maximumRHSError_;
+  CoinWorkDouble maximumRHSError_;
   /// maximumBoundInfeasibility.
-  double maximumBoundInfeasibility_;
+  CoinWorkDouble maximumBoundInfeasibility_;
   /// maximumDualError.
-  double maximumDualError_;
+  CoinWorkDouble maximumDualError_;
   /// diagonalScaleFactor.
-  double diagonalScaleFactor_;
+  CoinWorkDouble diagonalScaleFactor_;
   /// scaleFactor.  For scaling objective
-  double scaleFactor_;
+  CoinWorkDouble scaleFactor_;
   /// actualPrimalStep
-  double actualPrimalStep_;
+  CoinWorkDouble actualPrimalStep_;
   /// actualDualStep
-  double actualDualStep_;
+  CoinWorkDouble actualDualStep_;
   /// smallestInfeasibility
-  double smallestInfeasibility_;
+  CoinWorkDouble smallestInfeasibility_;
   /// historyInfeasibility.
 #define LENGTH_HISTORY 5
-  double historyInfeasibility_[LENGTH_HISTORY];
+  CoinWorkDouble historyInfeasibility_[LENGTH_HISTORY];
   /// complementarityGap.
-  double complementarityGap_;
+  CoinWorkDouble complementarityGap_;
   /// baseObjectiveNorm
-  double baseObjectiveNorm_;
+  CoinWorkDouble baseObjectiveNorm_;
   /// worstDirectionAccuracy
-  double worstDirectionAccuracy_;
+  CoinWorkDouble worstDirectionAccuracy_;
   /// maximumRHSChange
-  double maximumRHSChange_;
+  CoinWorkDouble maximumRHSChange_;
   /// errorRegion. i.e. Ax
-  double * errorRegion_;
+  CoinWorkDouble * errorRegion_;
   /// rhsFixRegion.
-  double * rhsFixRegion_;
+  CoinWorkDouble * rhsFixRegion_;
   /// upperSlack
-  double * upperSlack_;
+  CoinWorkDouble * upperSlack_;
   /// lowerSlack
-  double * lowerSlack_;
+  CoinWorkDouble * lowerSlack_;
   /// diagonal
-  double * diagonal_;
+  CoinWorkDouble * diagonal_;
   /// solution
-  double * solution_;
+  CoinWorkDouble * solution_;
   /// work array
-  double * workArray_;
+  CoinWorkDouble * workArray_;
   /// delta X
-  double * deltaX_;
+  CoinWorkDouble * deltaX_;
   /// delta Y
-  double * deltaY_;
+  CoinWorkDouble * deltaY_;
   /// deltaZ.
-  double * deltaZ_;
+  CoinWorkDouble * deltaZ_;
   /// deltaW.
-  double * deltaW_;
+  CoinWorkDouble * deltaW_;
   /// deltaS.
-  double * deltaSU_;
-  double * deltaSL_;
+  CoinWorkDouble * deltaSU_;
+  CoinWorkDouble * deltaSL_;
   /// Primal regularization array
-  double * primalR_;
+  CoinWorkDouble * primalR_;
   /// Dual regularization array
-  double * dualR_;
+  CoinWorkDouble * dualR_;
   /// rhs B
-  double * rhsB_;
+  CoinWorkDouble * rhsB_;
   /// rhsU.
-  double * rhsU_;
+  CoinWorkDouble * rhsU_;
   /// rhsL.
-  double * rhsL_;
+  CoinWorkDouble * rhsL_;
   /// rhsZ.
-  double * rhsZ_;
+  CoinWorkDouble * rhsZ_;
   /// rhsW.
-  double * rhsW_;
+  CoinWorkDouble * rhsW_;
   /// rhs C
-  double * rhsC_;
+  CoinWorkDouble * rhsC_;
   /// zVec
-  double * zVec_;
+  CoinWorkDouble * zVec_;
   /// wVec
-  double * wVec_;
+  CoinWorkDouble * wVec_;
   /// cholesky.
   ClpCholeskyBase * cholesky_;

trunk/Clp/src/ClpMatrixBase.cpp

@@ -605 +605 @@
   abort();
 }
+#if COIN_LONG_WORK 
+// For long double versions (aborts if not supported)
+void 
+ClpMatrixBase::times(CoinWorkDouble scalar,
+      const CoinWorkDouble * x, CoinWorkDouble * y) const 
+{ 
+  std::cerr<<"long times not supported - ClpMatrixBase"<<std::endl;
+  abort();
+}
+void 
+ClpMatrixBase::transposeTimes(CoinWorkDouble scalar,
+                              const CoinWorkDouble * x, CoinWorkDouble * y) const 
+{ 
+  std::cerr<<"long transposeTimes not supported - ClpMatrixBase"<<std::endl;
+  abort();
+}
+#endif

trunk/Clp/src/ClpMatrixBase.hpp

@@ -5 +5 @@
 
 #include "CoinPragma.hpp"
+#include "CoinFinite.hpp"
 
 #include "CoinPackedMatrix.hpp"
@@ -277 +278 @@
                               const double * columnScale,
                               double * spare=NULL) const;
+#if COIN_LONG_WORK 
+  // For long double versions (aborts if not supported)
+  virtual void times(CoinWorkDouble scalar,
+                     const CoinWorkDouble * x, CoinWorkDouble * y) const ;
+  virtual void transposeTimes(CoinWorkDouble scalar,
+                              const CoinWorkDouble * x, CoinWorkDouble * y) const ;
+#endif
   /** Return <code>x * scalar *A + y</code> in <code>z</code>. 
       Can use y as temporary array (will be empty at end)

trunk/Clp/src/ClpMessage.cpp

@@ -82 +82 @@
   {CLP_BARRIER_INFO,45,3,"Detail - %s"}, 
   {CLP_BARRIER_END,46,1,"At end primal/dual infeasibilities %g/%g, complementarity gap %g, objective %g"}, 
-  {CLP_BARRIER_ACCURACY,47,2,"Relative error in phase %d, refinement %d is %g"}, 
+  {CLP_BARRIER_ACCURACY,47,2,"Relative error in phase %d, %d refinements reduced from %g to %g"}, 
   {CLP_BARRIER_SAFE,48,2,"Initial safe primal value %g, objective norm %g"}, 
   {CLP_BARRIER_NEGATIVE_GAPS,49,3,"%d negative gaps summing to %g"}, 

trunk/Clp/src/ClpPackedMatrix.cpp

@@ -5961 +5961 @@
   output->setPackedMode(true);
 }
+#if COIN_LONG_WORK 
+// For long double versions
+void 
+ClpPackedMatrix::times(CoinWorkDouble scalar,
+                   const CoinWorkDouble * x, CoinWorkDouble * y) const
+{
+  int iRow,iColumn;
+  // get matrix data pointers
+  const int * row = matrix_->getIndices();
+  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
+  const double * elementByColumn = matrix_->getElements();
+  //memset(y,0,matrix_->getNumRows()*sizeof(double));
+  assert (((flags_&2)!=0)==matrix_->hasGaps());
+  if (!(flags_&2)) {
+    for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
+      CoinBigIndex j;
+      CoinWorkDouble value = x[iColumn];
+      if (value) {
+        CoinBigIndex start = columnStart[iColumn];
+        CoinBigIndex end = columnStart[iColumn+1];
+        value *= scalar;
+        for (j=start; j<end; j++) {
+          iRow=row[j];
+          y[iRow] += value*elementByColumn[j];
+        }
+      }
+    }
+  } else {
+    const int * columnLength = matrix_->getVectorLengths(); 
+    for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
+      CoinBigIndex j;
+      CoinWorkDouble value = x[iColumn];
+      if (value) {
+        CoinBigIndex start = columnStart[iColumn];
+        CoinBigIndex end = start + columnLength[iColumn];
+        value *= scalar;
+        for (j=start; j<end; j++) {
+          iRow=row[j];
+          y[iRow] += value*elementByColumn[j];
+        }
+      }
+    }
+  }
+}
+void 
+ClpPackedMatrix::transposeTimes(CoinWorkDouble scalar,
+                                const CoinWorkDouble * x, CoinWorkDouble * y) const
+{
+  int iColumn;
+  // get matrix data pointers
+  const int * row = matrix_->getIndices();
+  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
+  const double * elementByColumn = matrix_->getElements();
+  if (!(flags_&2)) {
+    if (scalar==-1.0) {
+      CoinBigIndex start=columnStart[0];
+      for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
+        CoinBigIndex j;
+        CoinBigIndex next=columnStart[iColumn+1];
+        CoinWorkDouble value=y[iColumn];
+        for (j=start;j<next;j++) {
+          int jRow=row[j];
+          value -= x[jRow]*elementByColumn[j];
+        }
+        start=next;
+        y[iColumn] = value;
+      }
+    } else {
+      CoinBigIndex start=columnStart[0];
+      for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
+        CoinBigIndex j;
+        CoinBigIndex next=columnStart[iColumn+1];
+        CoinWorkDouble value=0.0;
+        for (j=start;j<next;j++) {
+          int jRow=row[j];
+          value += x[jRow]*elementByColumn[j];
+        }
+        start=next;
+        y[iColumn] += value*scalar;
+      }
+    }
+  } else {
+    const int * columnLength = matrix_->getVectorLengths(); 
+    for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
+      CoinBigIndex j;
+      CoinWorkDouble value=0.0;
+      CoinBigIndex start = columnStart[iColumn];
+      CoinBigIndex end = start + columnLength[iColumn];
+      for (j=start; j<end; j++) {
+        int jRow=row[j];
+        value += x[jRow]*elementByColumn[j];
+      }
+      y[iColumn] += value*scalar;
+    }
+  }
+}
+#endif
 #ifdef CLP_ALL_ONE_FILE
 #undef reference

trunk/Clp/src/ClpPackedMatrix.hpp

@@ -264 +264 @@
   /// Sets up an effective RHS
   void useEffectiveRhs(ClpSimplex * model);
+#if COIN_LONG_WORK 
+  // For long double versions 
+  virtual void times(CoinWorkDouble scalar,
+                     const CoinWorkDouble * x, CoinWorkDouble * y) const ;
+  virtual void transposeTimes(CoinWorkDouble scalar,
+                              const CoinWorkDouble * x, CoinWorkDouble * y) const ;
+#endif
 //@}
 

trunk/Clp/src/ClpPredictorCorrector.cpp

@@ -33 +33 @@
     fwrite(&numberRows,sizeof(int),1,fp);
     fwrite(&numberColumns,sizeof(int),1,fp);
-    double dsave[20];
+    CoinWorkDouble dsave[20];
     memset(dsave,0,sizeof(dsave));
-    fwrite(dsave,sizeof(double),20,fp);
+    fwrite(dsave,sizeof(CoinWorkDouble),20,fp);
     int msave[20];
     memset(msave,0,sizeof(msave));
     msave[0]=numberIterations_;
     fwrite(msave,sizeof(int),20,fp);
-    fwrite(dual_,sizeof(double),numberRows,fp);
-    fwrite(errorRegion_,sizeof(double),numberRows,fp);
-    fwrite(rhsFixRegion_,sizeof(double),numberRows,fp);
-    fwrite(solution_,sizeof(double),numberColumns,fp);
-    fwrite(solution_+numberColumns,sizeof(double),numberRows,fp);
-    fwrite(diagonal_,sizeof(double),numberColumns,fp);
-    fwrite(diagonal_+numberColumns,sizeof(double),numberRows,fp);
-    fwrite(wVec_,sizeof(double),numberColumns,fp);
-    fwrite(wVec_+numberColumns,sizeof(double),numberRows,fp);
-    fwrite(zVec_,sizeof(double),numberColumns,fp);
-    fwrite(zVec_+numberColumns,sizeof(double),numberRows,fp);
-    fwrite(upperSlack_,sizeof(double),numberColumns,fp);
-    fwrite(upperSlack_+numberColumns,sizeof(double),numberRows,fp);
-    fwrite(lowerSlack_,sizeof(double),numberColumns,fp);
-    fwrite(lowerSlack_+numberColumns,sizeof(double),numberRows,fp);
+    fwrite(dual_,sizeof(CoinWorkDouble),numberRows,fp);
+    fwrite(errorRegion_,sizeof(CoinWorkDouble),numberRows,fp);
+    fwrite(rhsFixRegion_,sizeof(CoinWorkDouble),numberRows,fp);
+    fwrite(solution_,sizeof(CoinWorkDouble),numberColumns,fp);
+    fwrite(solution_+numberColumns,sizeof(CoinWorkDouble),numberRows,fp);
+    fwrite(diagonal_,sizeof(CoinWorkDouble),numberColumns,fp);
+    fwrite(diagonal_+numberColumns,sizeof(CoinWorkDouble),numberRows,fp);
+    fwrite(wVec_,sizeof(CoinWorkDouble),numberColumns,fp);
+    fwrite(wVec_+numberColumns,sizeof(CoinWorkDouble),numberRows,fp);
+    fwrite(zVec_,sizeof(CoinWorkDouble),numberColumns,fp);
+    fwrite(zVec_+numberColumns,sizeof(CoinWorkDouble),numberRows,fp);
+    fwrite(upperSlack_,sizeof(CoinWorkDouble),numberColumns,fp);
+    fwrite(upperSlack_+numberColumns,sizeof(CoinWorkDouble),numberRows,fp);
+    fwrite(lowerSlack_,sizeof(CoinWorkDouble),numberColumns,fp);
+    fwrite(lowerSlack_+numberColumns,sizeof(CoinWorkDouble),numberRows,fp);
     fclose(fp);
   } else {
@@ -61 +61 @@
 }
 #endif
-static double eScale=1.0e27;
-static double eBaseCaution=1.0e-12;
-static double eBase=1.0e-12;
-static double eRatio=1.0e40;
-static double eRatioCaution=1.0e25;
-static double eDiagonal=1.0e25;
-static double eDiagonalCaution=1.0e18;
-static double eExtra=1.0e-12;
+// Could change on CLP_LONG_CHOLESKY or COIN_LONG_WORK?
+static CoinWorkDouble eScale=1.0e27;
+static CoinWorkDouble eBaseCaution=1.0e-12;
+static CoinWorkDouble eBase=1.0e-12;
+static CoinWorkDouble eRatio=1.0e40;
+static CoinWorkDouble eRatioCaution=1.0e25;
+static CoinWorkDouble eDiagonal=1.0e25;
+static CoinWorkDouble eDiagonalCaution=1.0e18;
+static CoinWorkDouble eExtra=1.0e-12;
 
 // main function
@@ -81 +82 @@
     return 2;
   }
+#if COIN_LONG_WORK 
+  // reallocate some regions
+  double * dualSave = dual_;
+  dual_ = reinterpret_cast<double *>(new CoinWorkDouble[numberRows_]);
+  double * reducedCostSave = reducedCost_;
+  reducedCost_ = reinterpret_cast<double *>(new CoinWorkDouble[numberColumns_]);
+#endif  
   //diagonalPerturbation_=1.0e-25;
   ClpMatrixBase * saveMatrix = NULL;
@@ -165 +173 @@
     return -1;
   }
+  CoinWorkDouble * dualArray = reinterpret_cast<CoinWorkDouble *>(dual_);
   // Could try centering steps without any original step i.e. just center
   //firstFactorization(false);
-  CoinZeroN(dual_,numberRows_);
+  CoinZeroN(dualArray,numberRows_);
   multiplyAdd(solution_+numberColumns_,numberRows_,-1.0,errorRegion_,0.0);
   matrix_->times(1.0,solution_,errorRegion_);
   maximumRHSError_=maximumAbsElement(errorRegion_,numberRows_);
   maximumBoundInfeasibility_=maximumRHSError_;
-  //double maximumDualError_=COIN_DBL_MAX;
+  //CoinWorkDouble maximumDualError_=COIN_DBL_MAX;
   //initialize
   actualDualStep_=0.0;
@@ -184 +193 @@
   int numberFixedTotal=numberFixed;
   //int numberRows_DroppedBefore=0;
-  //double extra=eExtra;
-  //double maximumPerturbation=COIN_DBL_MAX;
+  //CoinWorkDouble extra=eExtra;
+  //CoinWorkDouble maximumPerturbation=COIN_DBL_MAX;
   //constants for infeas interior point
-  const double beta2 = 0.99995;
-  const double tau   = 0.00002;
-  double lastComplementarityGap=COIN_DBL_MAX;
-  double lastStep=1.0;
+  const CoinWorkDouble beta2 = 0.99995;
+  const CoinWorkDouble tau   = 0.00002;
+  CoinWorkDouble lastComplementarityGap=COIN_DBL_MAX;
+  CoinWorkDouble lastStep=1.0;
   // use to see if to take affine
-  double checkGap = COIN_DBL_MAX;
+  CoinWorkDouble checkGap = COIN_DBL_MAX;
   int lastGoodIteration=0;
-  double bestObjectiveGap=COIN_DBL_MAX;
-  double bestObjective=COIN_DBL_MAX;
+  CoinWorkDouble bestObjectiveGap=COIN_DBL_MAX;
+  CoinWorkDouble bestObjective=COIN_DBL_MAX;
+  int bestKilled=-1;
   int saveIteration=-1;
+  int saveIteration2=-1;
   bool sloppyOptimal=false;
-  double * savePi=NULL;
-  double * savePrimal=NULL;
+  CoinWorkDouble * savePi=NULL;
+  CoinWorkDouble * savePrimal=NULL;
+  CoinWorkDouble * savePi2=NULL;
+  CoinWorkDouble * savePrimal2=NULL;
   // Extra regions for centering
-  double * saveX = new double[numberTotal];
-  double * saveY = new double[numberRows_];
-  double * saveZ = new double[numberTotal];
-  double * saveW = new double[numberTotal];
-  double * saveSL = new double[numberTotal];
-  double * saveSU = new double[numberTotal];
+  CoinWorkDouble * saveX = new CoinWorkDouble[numberTotal];
+  CoinWorkDouble * saveY = new CoinWorkDouble[numberRows_];
+  CoinWorkDouble * saveZ = new CoinWorkDouble[numberTotal];
+  CoinWorkDouble * saveW = new CoinWorkDouble[numberTotal];
+  CoinWorkDouble * saveSL = new CoinWorkDouble[numberTotal];
+  CoinWorkDouble * saveSU = new CoinWorkDouble[numberTotal];
   // Save smallest mu used in primal dual moves
-  double smallestPrimalDualMu=COIN_DBL_MAX;
-  double objScale = optimizationDirection_/
+  CoinWorkDouble smallestPrimalDualMu=COIN_DBL_MAX;
+  CoinWorkDouble objScale = optimizationDirection_/
     (rhsScale_*objectiveScale_);
   while (problemStatus_<0) {
@@ -231 +244 @@
     handler_->message(CLP_BARRIER_ITERATION,messages_)
       <<numberIterations_
-      <<primalObjective_*objScale- dblParam_[ClpObjOffset]
-      << dualObjective_*objScale- dblParam_[ClpObjOffset]
-      <<complementarityGap_
+      <<static_cast<double>(primalObjective_*objScale- dblParam_[ClpObjOffset])
+      << static_cast<double>(dualObjective_*objScale- dblParam_[ClpObjOffset])
+      <<static_cast<double>(complementarityGap_)
       <<numberFixedTotal
       <<cholesky_->rank()
@@ -253 +266 @@
     //saveIteration=-1;
     lastStep = CoinMin(actualPrimalStep_,actualDualStep_);
-    double goodGapChange;
-    //#define KEEP_GOING_IF_FIXED 0
+    CoinWorkDouble goodGapChange;
+    //#define KEEP_GOING_IF_FIXED 5
 #ifndef KEEP_GOING_IF_FIXED 
 #define KEEP_GOING_IF_FIXED 10000
@@ -265 +278 @@
         goodGapChange=0.99; // make more likely to carry on
     } 
-    double gapO;
-    double lastGood=bestObjectiveGap;
+    CoinWorkDouble gapO;
+    CoinWorkDouble lastGood=bestObjectiveGap;
     if (gonePrimalFeasible_&&goneDualFeasible_) {
-      double largestObjective;
-      if (fabs(primalObjective_)>fabs(dualObjective_)) {
-        largestObjective = fabs(primalObjective_);
+      CoinWorkDouble largestObjective;
+      if (CoinAbs(primalObjective_)>CoinAbs(dualObjective_)) {
+        largestObjective = CoinAbs(primalObjective_);
       } else {
-        largestObjective = fabs(dualObjective_);
+        largestObjective = CoinAbs(dualObjective_);
       } 
       if (largestObjective<1.0) {
         largestObjective=1.0;
       } 
-      gapO=fabs(primalObjective_-dualObjective_)/largestObjective;
+      gapO=CoinAbs(primalObjective_-dualObjective_)/largestObjective;
       handler_->message(CLP_BARRIER_OBJECTIVE_GAP,messages_)
-        <<gapO
+        <<static_cast<double>(gapO)
         <<CoinMessageEol;
       //start saving best
-      if (gapO<bestObjectiveGap) 
+      bool saveIt=false;
+      if (gapO<bestObjectiveGap) { 
         bestObjectiveGap=gapO;
+#ifndef SAVE_ON_OBJ
+        saveIt=true;
+#endif
+      }
       if (primalObjective_<bestObjective) {
+        bestObjective=primalObjective_;
+#ifdef SAVE_ON_OBJ
+        saveIt=true;
+#endif
+      }
+      if (numberFixedTotal>bestKilled) {
+        bestKilled=numberFixedTotal;
+#if KEEP_GOING_IF_FIXED<10
+        saveIt=true;
+#endif
+      }
+      if (saveIt) {
+#if KEEP_GOING_IF_FIXED<10
+        printf("saving\n");
+#endif
         saveIteration=numberIterations_;
-        bestObjective=primalObjective_;
         if (!savePi) {
-          savePi=new double[numberRows_];
-          savePrimal = new double [numberTotal];
-        } 
-        CoinMemcpyN(dual_,numberRows_,savePi);
+          savePi=new CoinWorkDouble[numberRows_];
+          savePrimal = new CoinWorkDouble [numberTotal];
+        } 
+        CoinMemcpyN(dualArray,numberRows_,savePi);
         CoinMemcpyN(solution_,numberTotal,savePrimal);
       } else if(gapO>2.0*bestObjectiveGap) {
@@ -300 +332 @@
       //std::cout <<"could stop"<<std::endl;
       //gapO=0.0;
-      if (fabs(primalObjective_-dualObjective_)<dualTolerance()) {
+      if (CoinAbs(primalObjective_-dualObjective_)<dualTolerance()) {
         gapO=0.0;
       } 
@@ -308 +340 @@
         handler_->message(CLP_BARRIER_GONE_INFEASIBLE,messages_)
           <<CoinMessageEol;
+        CoinWorkDouble scaledRHSError=maximumRHSError_/(solutionNorm_+10.0);
+        // save alternate
+        if (numberFixedTotal>bestKilled&&
+            maximumBoundInfeasibility_<1.0e-6&&
+            scaledRHSError<1.0e-2) {
+          bestKilled=numberFixedTotal;
+#if KEEP_GOING_IF_FIXED<10
+          printf("saving alternate\n");
+#endif
+          saveIteration2=numberIterations_;
+          if (!savePi2) {
+            savePi2=new CoinWorkDouble[numberRows_];
+            savePrimal2 = new CoinWorkDouble [numberTotal];
+          } 
+          CoinMemcpyN(dualArray,numberRows_,savePi2);
+          CoinMemcpyN(solution_,numberTotal,savePrimal2);
+        }
         if (sloppyOptimal) {
           // vaguely optimal
-          double scaledRHSError=maximumRHSError_/solutionNorm_;
           if (maximumBoundInfeasibility_>1.0e-2||
               scaledRHSError>1.0e-2||
@@ -322 +370 @@
         } else {
           // not close to optimal but check if getting bad
-          double scaledRHSError=maximumRHSError_/solutionNorm_;
+          CoinWorkDouble scaledRHSError=maximumRHSError_/(solutionNorm_+10.0);
           if ((maximumBoundInfeasibility_>1.0e-1||
               scaledRHSError>1.0e-1||
@@ -354 +402 @@
       sloppyOptimal=true;
       handler_->message(CLP_BARRIER_CLOSE_TO_OPTIMAL,messages_)
-        <<numberIterations_<<complementarityGap_
+        <<numberIterations_<<static_cast<double>(complementarityGap_)
         <<CoinMessageEol;
     } 
@@ -363 +411 @@
     if (complementarityGap_>=1.05*lastComplementarityGap) {
       handler_->message(CLP_BARRIER_COMPLEMENTARITY,messages_)
-        <<complementarityGap_<<"increasing"
+        <<static_cast<double>(complementarityGap_)<<"increasing"
         <<CoinMessageEol;
       if (saveIteration>=0&&sloppyOptimal) {
@@ -380 +428 @@
                complementarityGap_<1.0e-3) {
       handler_->message(CLP_BARRIER_COMPLEMENTARITY,messages_)
-        <<complementarityGap_<<"not decreasing"
+        <<static_cast<double>(complementarityGap_)<<"not decreasing"
         <<CoinMessageEol;
       if (gapO>0.75*lastGood&&numberFixed<KEEP_GOING_IF_FIXED) {
@@ -388 +436 @@
                complementarityGap_<1.0e-6) {
       handler_->message(CLP_BARRIER_COMPLEMENTARITY,messages_)
-        <<complementarityGap_<<"not decreasing"
+        <<static_cast<double>(complementarityGap_)<<"not decreasing"
         <<CoinMessageEol;
       break;
@@ -421 +469 @@
     } 
     if (gapO<1.0e-9) {
-      double value=gapO*complementarityGap_;
+      CoinWorkDouble value=gapO*complementarityGap_;
       value*=actualPrimalStep_;
       value*=actualDualStep_;
@@ -433 +481 @@
       } 
     } 
-    double nextGap=COIN_DBL_MAX;
+    CoinWorkDouble nextGap=COIN_DBL_MAX;
     int nextNumber=0;
     int nextNumberItems=0;
@@ -441 +489 @@
     //prepare for cholesky.  Set up scaled diagonal in deltaX
     //  ** for efficiency may be better if scale factor known before
-    double norm2=0.0;
-    double maximumValue;
+    CoinWorkDouble norm2=0.0;
+    CoinWorkDouble maximumValue;
     getNorms(diagonal_,numberTotal,maximumValue,norm2);
-    diagonalNorm_ = sqrt(norm2/numberComplementarityPairs_);
+    diagonalNorm_ = CoinSqrt(norm2/numberComplementarityPairs_);
     diagonalScaleFactor_=1.0;
-    double maximumAllowable=eScale;
+    CoinWorkDouble maximumAllowable=eScale;
     //scale so largest is less than allowable ? could do better
-    double factor=0.5;
+    CoinWorkDouble factor=0.5;
     while (maximumValue>maximumAllowable) {
       diagonalScaleFactor_*=factor;
@@ -455 +503 @@
     if (diagonalScaleFactor_!=1.0) {
       handler_->message(CLP_BARRIER_SCALING,messages_)
-        <<"diagonal"<<diagonalScaleFactor_
+        <<"diagonal"<<static_cast<double>(diagonalScaleFactor_)
         <<CoinMessageEol;
       diagonalNorm_*=diagonalScaleFactor_;
@@ -496 +544 @@
     }
     int phase=0; // predictor, corrector , primal dual
-    double directionAccuracy=0.0;
+    CoinWorkDouble directionAccuracy=0.0;
     bool doCorrector=true;
     bool goodMove=true;
@@ -517 +565 @@
       debugMove(0,1.0e-2,1.0e-2);
     }
-    double affineGap=nextGap;
+    CoinWorkDouble affineGap=nextGap;
     int bestPhase=0;
-    double bestNextGap=nextGap;
+    CoinWorkDouble bestNextGap=nextGap;
     // ?
     bestNextGap=CoinMax(nextGap,0.8*complementarityGap_);
@@ -526 +574 @@
     if (complementarityGap_>1.0e-4*numberComplementarityPairs_) {
       //std::cout <<"predicted duality gap "<<nextGap<<std::endl;
-      double part1=nextGap/numberComplementarityPairs_;
+      CoinWorkDouble part1=nextGap/numberComplementarityPairs_;
       part1=nextGap/numberComplementarityItems_;
-      double part2=nextGap/complementarityGap_;
+      CoinWorkDouble part2=nextGap/complementarityGap_;
       mu_=part1*part2*part2;
 #if 0
-      double papermu =complementarityGap_/numberComplementarityPairs_;
-      double affmu = nextGap/nextNumber;
-      double sigma = pow(affmu/papermu,3);
+      CoinWorkDouble papermu =complementarityGap_/numberComplementarityPairs_;
+      CoinWorkDouble affmu = nextGap/nextNumber;
+      CoinWorkDouble sigma = pow(affmu/papermu,3);
       printf("mu %g, papermu %g, affmu %g, sigma %g sigmamu %g\n",
              mu_,papermu,affmu,sigma,sigma*papermu);
 #endif
       //printf("paper mu %g\n",(nextGap*nextGap*nextGap)/(complementarityGap_*complementarityGap_*
-      //                                            (double) numberComplementarityPairs_));
+      //                                            (CoinWorkDouble) numberComplementarityPairs_));
     } else {
-      double phi;
+      CoinWorkDouble phi;
       if (numberComplementarityPairs_<=5000) {
-        phi=pow(static_cast<double> (numberComplementarityPairs_),2.0);
+        phi=pow(static_cast<CoinWorkDouble> (numberComplementarityPairs_),2.0);
       } else {
-        phi=pow(static_cast<double> (numberComplementarityPairs_),1.5);
+        phi=pow(static_cast<CoinWorkDouble> (numberComplementarityPairs_),1.5);
         if (phi<500.0*500.0) {
           phi=500.0*500.0;
@@ -552 +600 @@
     } 
     //save information
-    double product=affineProduct();
+    CoinWorkDouble product=affineProduct();
 #if 0
     //can we do corrector step?
-    double xx= complementarityGap_*(beta2-tau) +product;
+    CoinWorkDouble xx= complementarityGap_*(beta2-tau) +product;
     if (xx>0.0) {
-      double saveMu = mu_;
-      double mu2=numberComplementarityPairs_;
+      CoinWorkDouble saveMu = mu_;
+      CoinWorkDouble mu2=numberComplementarityPairs_;
       mu2=xx/mu2;
       if (mu2>mu_) {
@@ -593 +641 @@
         bestNextGap=COIN_DBL_MAX;
       }
-      //double floatNumber = 2.0*numberComplementarityPairs_;
+      //CoinWorkDouble floatNumber = 2.0*numberComplementarityPairs_;
       //floatNumber = 1.0*numberComplementarityItems_;
       //mu_=nextGap/floatNumber;
@@ -615 +663 @@
       CoinMemcpyN(deltaSU_,numberTotal,saveSU);
 #ifdef HALVE
-      double savePrimalStep = actualPrimalStep_;
-      double saveDualStep = actualDualStep_;
-      double saveMu = mu_;
+      CoinWorkDouble savePrimalStep = actualPrimalStep_;
+      CoinWorkDouble saveDualStep = actualDualStep_;
+      CoinWorkDouble saveMu = mu_;
 #endif
       //set up for next step
       setupForSolve(phase);
-      double directionAccuracy2=findDirectionVector(phase);
+      CoinWorkDouble directionAccuracy2=findDirectionVector(phase);
       if (directionAccuracy2>worstDirectionAccuracy_) {
         worstDirectionAccuracy_=directionAccuracy2;
       } 
-      double testValue=1.0e2*directionAccuracy;
+      CoinWorkDouble testValue=1.0e2*directionAccuracy;
       if (1.0e2*projectionTolerance_>testValue) {
         testValue=1.0e2*projectionTolerance_;
@@ -644 +692 @@
       if (goodMove) {
         phase=1;
-        double norm = findStepLength(phase);
+        CoinWorkDouble norm = findStepLength(phase);
         nextGap = complementarityGap(nextNumber,nextNumberItems,1);
         debugMove(1,actualPrimalStep_,actualDualStep_);
@@ -668 +716 @@
         //printf("halve %d\n",nHalve);
         nHalve++;
-        const double lambda=0.5;
+        const CoinWorkDouble lambda=0.5;
         for (i=0;i<numberRows_;i++)
           deltaY_[i] = lambda*deltaY_[i]+(1.0-lambda)*saveY[i];
@@ -699 +747 @@
     if (!goodMove) {
       // Just primal dual step
-      double floatNumber;
+      CoinWorkDouble floatNumber;
       floatNumber = 2.0*numberComplementarityPairs_;
       //floatNumber = numberComplementarityItems_;
-      double saveMu=mu_; // use one from predictor corrector
+      CoinWorkDouble saveMu=mu_; // use one from predictor corrector
       mu_=complementarityGap_/floatNumber;
       // If going well try small mu
-      mu_ *= sqrt((1.0-lastStep)/(1.0+10.0*lastStep));
-      double mu1=mu_;
-      double phi;
+      mu_ *= CoinSqrt((1.0-lastStep)/(1.0+10.0*lastStep));
+      CoinWorkDouble mu1=mu_;
+      CoinWorkDouble phi;
       if (numberComplementarityPairs_<=500) {
-        phi=pow(static_cast<double> (numberComplementarityPairs_),2.0);
+        phi=pow(static_cast<CoinWorkDouble> (numberComplementarityPairs_),2.0);
       } else {
-        phi=pow(static_cast<double> (numberComplementarityPairs_),1.5);
+        phi=pow(static_cast<CoinWorkDouble> (numberComplementarityPairs_),1.5);
         if (phi<500.0*500.0) {
           phi=500.0*500.0;
@@ -719 +767 @@
       //printf("pd mu %g, alternate %g, smallest %g\n",
       //     mu_,mu1,smallestPrimalDualMu);
-      mu_ = sqrt(mu_*mu1);
+      mu_ = CoinSqrt(mu_*mu1);
       mu_=mu1;
       if ((numberIterations_&1)==0||numberIterations_<10)
@@ -737 +785 @@
       setupForSolve(2);
       findDirectionVector(2);
-      double norm=findStepLength(2);
+      CoinWorkDouble norm=findStepLength(2);
       // just for debug
       bestNextGap = complementarityGap_*1.0005;
@@ -760 +808 @@
       smallestPrimalDualMu=mu_;
     if (!goodMove)
-      mu_=nextGap / (static_cast<double> (nextNumber)*1.1);
+      mu_=nextGap / (static_cast<CoinWorkDouble> (nextNumber)*1.1);
     //if (quadraticObj)
     //goodMove=true; 
@@ -766 +814 @@
     // Do centering steps
     int numberTries=0;
-    double nextCenterGap=0.0;
+    CoinWorkDouble nextCenterGap=0.0;
     int numberGoodTries=0;
-    double originalDualStep=actualDualStep_;
-    double originalPrimalStep=actualPrimalStep_;
+    CoinWorkDouble originalDualStep=actualDualStep_;
+    CoinWorkDouble originalPrimalStep=actualPrimalStep_;
     if (actualDualStep_>0.9&&actualPrimalStep_>0.9)
       goodMove=false; // don't bother
@@ -781 +829 @@
       CoinMemcpyN(deltaZ_,numberTotal,saveZ);
       CoinMemcpyN(deltaW_,numberTotal,saveW);
-      double savePrimalStep = actualPrimalStep_;
-      double saveDualStep = actualDualStep_;
-      double saveMu = mu_;
+      CoinWorkDouble savePrimalStep = actualPrimalStep_;
+      CoinWorkDouble saveDualStep = actualDualStep_;
+      CoinWorkDouble saveMu = mu_;
       setupForSolve(3);
       findDirectionVector(3);
@@ -789 +837 @@
       debugMove(3,actualPrimalStep_,actualDualStep_);
       //debugMove(3,1.0e-7,1.0e-7);
-      double xGap = complementarityGap(nextNumber,nextNumberItems,3);
+      CoinWorkDouble xGap = complementarityGap(nextNumber,nextNumberItems,3);
       // If one small then that's the one that counts
-      double checkDual=saveDualStep;
-      double checkPrimal=savePrimalStep;
+      CoinWorkDouble checkDual=saveDualStep;
+      CoinWorkDouble checkPrimal=savePrimalStep;
       if (checkDual>5.0*checkPrimal) {
         checkDual=2.0*checkPrimal;
@@ -848 +896 @@
   delete [] saveSU;
   if (savePi) {
-    std::cout<<"Restoring from iteration "<<saveIteration<<std::endl;
-    CoinMemcpyN(savePi,numberRows_,dual_);
-    CoinMemcpyN(savePrimal,numberTotal,solution_);
+    if (numberIterations_-saveIteration>20&&
+        numberIterations_-saveIteration2<5) { 
+#if KEEP_GOING_IF_FIXED<10
+      std::cout<<"Restoring2 from iteration "<<saveIteration2<<std::endl;
+#endif
+      CoinMemcpyN(savePi2,numberRows_,dualArray);
+      CoinMemcpyN(savePrimal2,numberTotal,solution_);
+    } else {
+#if KEEP_GOING_IF_FIXED<10
+      std::cout<<"Restoring from iteration "<<saveIteration<<std::endl;
+#endif
+      CoinMemcpyN(savePi,numberRows_,dualArray);
+      CoinMemcpyN(savePrimal,numberTotal,solution_);
+    }
     delete [] savePi;
     delete [] savePrimal;
   } 
+  delete [] savePi2;
+  delete [] savePrimal2;
   //recompute slacks
   // Split out solution
@@ -861 +922 @@
   //unscale objective
   multiplyAdd(NULL,numberTotal,0.0,cost_,scaleFactor_);
-  multiplyAdd(NULL,numberRows_,0,dual_,scaleFactor_);
+  multiplyAdd(NULL,numberRows_,0,dualArray,scaleFactor_);
   checkSolution();
-  CoinMemcpyN(reducedCost_,numberColumns_,dj_);
+  //CoinMemcpyN(reducedCost_,numberColumns_,dj_);
   // If quadratic use last solution 
   // Restore quadratic objective if necessary
@@ -872 +933 @@
   }
   handler_->message(CLP_BARRIER_END,messages_)
-    <<sumPrimalInfeasibilities_
-    <<sumDualInfeasibilities_
-    <<complementarityGap_
-    <<objectiveValue()
+    <<static_cast<double>(sumPrimalInfeasibilities_)
+    <<static_cast<double>(sumDualInfeasibilities_)
+    <<static_cast<double>(complementarityGap_)
+    <<static_cast<double>(objectiveValue())
     <<CoinMessageEol;
   //#ifdef SOME_DEBUG
@@ -886 +947 @@
   //std::cout<<"Primal objective "<<objectiveValue()<<std::endl;
   //std::cout<<"maximum complementarity "<<worstComplementarity_<<std::endl;
+#if COIN_LONG_WORK 
+  // put back dual
+  delete [] dual_;
+  delete [] reducedCost_;
+  dual_ = dualSave;
+  reducedCost_=reducedCostSave;
+#endif
   //delete all temporary regions
   deleteWorkingData();
+#if KEEP_GOING_IF_FIXED<10
+#ifndef NDEBUG
+  {
+    static int kk=0;
+    char name[20];
+    sprintf(name,"save.sol.%d",kk);
+    kk++;
+    printf("saving to file %s\n",name);
+    FILE * fp = fopen(name,"wb");
+    int numberWritten;
+    numberWritten = fwrite(&numberColumns_,sizeof(int),1,fp);
+    assert (numberWritten==1);
+    numberWritten = fwrite(columnActivity_,sizeof(double),numberColumns_,fp);
+    assert (numberWritten==numberColumns_);
+    fclose(fp);
+  }
+#endif
+#endif
   if (saveMatrix) {
     // restore normal copy
@@ -899 +985 @@
 //         1 corrector
 //         2 primal dual
-double ClpPredictorCorrector::findStepLength( int phase)
+CoinWorkDouble ClpPredictorCorrector::findStepLength( int phase)
 {
-  double directionNorm=0.0;
-  double maximumPrimalStep=COIN_DBL_MAX;
-  double maximumDualStep=COIN_DBL_MAX;
+  CoinWorkDouble directionNorm=0.0;
+  CoinWorkDouble maximumPrimalStep=COIN_DBL_MAX;
+  CoinWorkDouble maximumDualStep=COIN_DBL_MAX;
   int numberTotal = numberRows_+numberColumns_;
-  double tolerance = 1.0e-12;
+#if CLP_LONG_CHOLESKY<2
+  CoinWorkDouble tolerance = 1.0e-12;
+#else
+  CoinWorkDouble tolerance = 1.0e-14;
+#endif
   int chosenPrimalSequence=-1;
   int chosenDualSequence=-1;
@@ -911 +1001 @@
   bool lowDual=false;
   // If done many iterations then allow to hit boundary
-  double hitTolerance;
+  CoinWorkDouble hitTolerance;
   //printf("objective norm %g\n",objectiveNorm_);
   if (numberIterations_<80||!gonePrimalFeasible_)
@@ -922 +1012 @@
   for (iColumn=0;iColumn<numberTotal;iColumn++) {
     if (!flagged(iColumn)) {
-      double directionElement=deltaX_[iColumn];
-      if (directionNorm<fabs(directionElement)) {
-        directionNorm=fabs(directionElement);
-      }
-      if (0) 
-      printf("%d %g %g %g %g %g %g %g %g %g %g %g\n",
-             iColumn,solution_[iColumn],deltaX_[iColumn],
-             lowerSlack_[iColumn],deltaSL_[iColumn],
-             upperSlack_[iColumn],deltaSU_[iColumn],
-             dj_[iColumn],
-             zVec_[iColumn],deltaZ_[iColumn],
-             wVec_[iColumn],deltaW_[iColumn]);
+      CoinWorkDouble directionElement=deltaX_[iColumn];
+      if (directionNorm<CoinAbs(directionElement)) {
+        directionNorm=CoinAbs(directionElement);
+      }
       if (lowerBound(iColumn)) {
-        double delta = - deltaSL_[iColumn];
-        double z1 = deltaZ_[iColumn];
-        double newZ = zVec_[iColumn]+z1;
+        CoinWorkDouble delta = - deltaSL_[iColumn];
+        CoinWorkDouble z1 = deltaZ_[iColumn];
+        CoinWorkDouble newZ = zVec_[iColumn]+z1;
         if (zVec_[iColumn]>tolerance) {
           if (zVec_[iColumn]<-z1*maximumDualStep) {
@@ -946 +1028 @@
         } 
         if (lowerSlack_[iColumn]<maximumPrimalStep*delta) {
-          double newStep=lowerSlack_[iColumn]/delta;
+          CoinWorkDouble newStep=lowerSlack_[iColumn]/delta;
           if (newStep>0.2||newZ<hitTolerance||delta>1.0e3||delta<=1.0e-6||dj_[iColumn]<hitTolerance) {
             maximumPrimalStep = newStep;
@@ -957 +1039 @@
       }
       if (upperBound(iColumn)) {
-        double delta = - deltaSU_[iColumn];;
-        double w1 = deltaW_[iColumn];
-        double newT = wVec_[iColumn]+w1;
+        CoinWorkDouble delta = - deltaSU_[iColumn];;
+        CoinWorkDouble w1 = deltaW_[iColumn];
+        CoinWorkDouble newT = wVec_[iColumn]+w1;
         if (wVec_[iColumn]>tolerance) {
           if (wVec_[iColumn]<-w1*maximumDualStep) {
@@ -968 +1050 @@
         } 
         if (upperSlack_[iColumn]<maximumPrimalStep*delta) {
-          double newStep=upperSlack_[iColumn]/delta;
+          CoinWorkDouble newStep=upperSlack_[iColumn]/delta;
           if (newStep>0.2||newT<hitTolerance||delta>1.0e3||delta<=1.0e-6||dj_[iColumn]>-hitTolerance) {
             maximumPrimalStep = newStep;
@@ -1026 +1108 @@
   if (quadraticObj) {
     // Use smaller unless very small
-    double smallerStep=CoinMin(actualDualStep_,actualPrimalStep_);
+    CoinWorkDouble smallerStep=CoinMin(actualDualStep_,actualPrimalStep_);
     if (smallerStep>0.0001) {
       actualDualStep_=smallerStep;
@@ -1042 +1124 @@
       // at first - just check better - if not 
       // Complementarity gap will be a*change*change + b*change +c
-      double a=0.0;
-      double b=0.0;
-      double c=0.0;
+      CoinWorkDouble a=0.0;
+      CoinWorkDouble b=0.0;
+      CoinWorkDouble c=0.0;
       /* SQUARE of dual infeasibility will be:
          square of dj - ......
       */
-      double aq=0.0;
-      double bq=0.0;
-      double cq=0.0;
-      double gamma2 = gamma_*gamma_; // gamma*gamma will be added to diagonal
-      double * linearDjChange = new double[numberTotal];
+      CoinWorkDouble aq=0.0;
+      CoinWorkDouble bq=0.0;
+      CoinWorkDouble cq=0.0;
+      CoinWorkDouble gamma2 = gamma_*gamma_; // gamma*gamma will be added to diagonal
+      CoinWorkDouble * linearDjChange = new CoinWorkDouble[numberTotal];
       CoinZeroN(linearDjChange,numberColumns_);
       multiplyAdd(deltaY_,numberRows_,1.0,linearDjChange+numberColumns_,0.0);
@@ -1060 +1142 @@
       const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
       const int * columnQuadraticLength = quadratic->getVectorLengths();
-      double * quadraticElement = quadratic->getMutableElements();
+      CoinWorkDouble * quadraticElement = quadratic->getMutableElements();
       for (iColumn=0;iColumn<numberTotal;iColumn++) {
-        double oldPrimal = solution_[iColumn];
+        CoinWorkDouble oldPrimal = solution_[iColumn];
         if (!flagged(iColumn)) {
           if (lowerBound(iColumn)) {
-            double change =oldPrimal+deltaX_[iColumn]-lowerSlack_[iColumn]-lower_[iColumn];
+            CoinWorkDouble change =oldPrimal+deltaX_[iColumn]-lowerSlack_[iColumn]-lower_[iColumn];
             c += lowerSlack_[iColumn]*zVec_[iColumn];
             b += lowerSlack_[iColumn]*deltaZ_[iColumn]+zVec_[iColumn]*change;
@@ -1071 +1153 @@
           }
           if (upperBound(iColumn)) {
-            double change =upper_[iColumn]-oldPrimal-deltaX_[iColumn]-upperSlack_[iColumn];
+            CoinWorkDouble change =upper_[iColumn]-oldPrimal-deltaX_[iColumn]-upperSlack_[iColumn];
             c += upperSlack_[iColumn]*wVec_[iColumn];
             b += upperSlack_[iColumn]*deltaW_[iColumn]+wVec_[iColumn]*change;
@@ -1077 +1159 @@
           } 
           // new djs are dj_ + change*value
-          double djChange = linearDjChange[iColumn];
+          CoinWorkDouble djChange = linearDjChange[iColumn];
           if (iColumn<numberColumns_) {
             for (CoinBigIndex j=columnQuadraticStart[iColumn];
                  j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
               int jColumn = columnQuadratic[j];
-              double changeJ = deltaX_[jColumn];
-              double elementValue = quadraticElement[j];
+              CoinWorkDouble changeJ = deltaX_[jColumn];
+              CoinWorkDouble elementValue = quadraticElement[j];
               djChange += changeJ*elementValue;
             }
           }
-          double gammaTerm = gamma2;
+          CoinWorkDouble gammaTerm = gamma2;
           if (primalR_) {
             gammaTerm += primalR_[iColumn];
@@ -1093 +1175 @@
           djChange += gammaTerm;
           // and dual infeasibility
-          double oldInf = dj_[iColumn]-zVec_[iColumn]+wVec_[iColumn]+
+          CoinWorkDouble oldInf = dj_[iColumn]-zVec_[iColumn]+wVec_[iColumn]+
             gammaTerm*solution_[iColumn];
-          double changeInf = djChange-deltaZ_[iColumn]+deltaW_[iColumn];
+          CoinWorkDouble changeInf = djChange-deltaZ_[iColumn]+deltaW_[iColumn];
           cq += oldInf*oldInf;
           bq += 2.0*oldInf*changeInf;
@@ -1108 +1190 @@
           } 
           // new djs are dj_ + change*value
-          double djChange = linearDjChange[iColumn];
+          CoinWorkDouble djChange = linearDjChange[iColumn];
           if (iColumn<numberColumns_) {
             for (CoinBigIndex j=columnQuadraticStart[iColumn];
                  j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
               int jColumn = columnQuadratic[j];
-              double changeJ = deltaX_[jColumn];
-              double elementValue = quadraticElement[j];
+              CoinWorkDouble changeJ = deltaX_[jColumn];
+              CoinWorkDouble elementValue = quadraticElement[j];
               djChange += changeJ*elementValue;
             }
           }
-          double gammaTerm = gamma2;
+          CoinWorkDouble gammaTerm = gamma2;
           if (primalR_) {
             gammaTerm += primalR_[iColumn];
@@ -1124 +1206 @@
           djChange += gammaTerm;
           // and dual infeasibility
-          double oldInf = dj_[iColumn]-zVec_[iColumn]+wVec_[iColumn]+
+          CoinWorkDouble oldInf = dj_[iColumn]-zVec_[iColumn]+wVec_[iColumn]+
             gammaTerm*solution_[iColumn];
-          double changeInf = djChange-deltaZ_[iColumn]+deltaW_[iColumn];
+          CoinWorkDouble changeInf = djChange-deltaZ_[iColumn]+deltaW_[iColumn];
           cq += oldInf*oldInf;
           bq += 2.0*oldInf*changeInf;
@@ -1136 +1218 @@
       // maybe use inf and do line search
       // To check see if matches at current step
-      double step=actualPrimalStep_;
-      //Current gap + solutionNorm_ * sqrt (sum square inf)
-      double multiplier = solutionNorm_;
+      CoinWorkDouble step=actualPrimalStep_;
+      //Current gap + solutionNorm_ * CoinSqrt (sum square inf)
+      CoinWorkDouble multiplier = solutionNorm_;
       multiplier *= 0.01;
       multiplier=1.0;
-      double currentInf =  multiplier*sqrt(cq);
-      double nextInf =  multiplier*sqrt(CoinMax(cq+step*bq+step*step*aq,0.0));
-      double allowedIncrease=1.4;
+      CoinWorkDouble currentInf =  multiplier*CoinSqrt(cq);
+      CoinWorkDouble nextInf =  multiplier*CoinSqrt(CoinMax(cq+step*bq+step*step*aq,0.0));
+      CoinWorkDouble allowedIncrease=1.4;
 #ifdef SOME_DEBUG
       printf("lin %g %g %g -> %g\n",a,b,c,
@@ -1161 +1243 @@
         //cq = CoinMax(cq,10.0);
         // convert to (x+q)*(x+q) = w
-        double q = bq/(1.0*aq);
-        double w = CoinMax(q*q + (cq/aq)*(allowedIncrease-1.0),0.0);
-        w = sqrt(w);
-        double stepX = w-q;
+        CoinWorkDouble q = bq/(1.0*aq);
+        CoinWorkDouble w = CoinMax(q*q + (cq/aq)*(allowedIncrease-1.0),0.0);
+        w = CoinSqrt(w);
+        CoinWorkDouble stepX = w-q;
         step=stepX;
         nextInf = 
-          multiplier*sqrt(CoinMax(cq+step*bq+step*step*aq,0.0));
+          multiplier*CoinSqrt(CoinMax(cq+step*bq+step*step*aq,0.0));
 #ifdef SOME_DEBUG
         printf ("with step of %g dualInf is %g\n",
@@ -1180 +1262 @@
     // minimize complementarity 
     // Complementarity gap will be a*change*change + b*change +c
-    double a=0.0;
-    double b=0.0;
-    double c=0.0;
+    CoinWorkDouble a=0.0;
+    CoinWorkDouble b=0.0;
+    CoinWorkDouble c=0.0;
     for (iColumn=0;iColumn<numberTotal;iColumn++) {
-      double oldPrimal = solution_[iColumn];
+      CoinWorkDouble oldPrimal = solution_[iColumn];
       if (!flagged(iColumn)) {
         if (lowerBound(iColumn)) {
-          double change =oldPrimal+deltaX_[iColumn]-lowerSlack_[iColumn]-lower_[iColumn];
+          CoinWorkDouble change =oldPrimal+deltaX_[iColumn]-lowerSlack_[iColumn]-lower_[iColumn];
           c += lowerSlack_[iColumn]*zVec_[iColumn];
           b += lowerSlack_[iColumn]*deltaZ_[iColumn]+zVec_[iColumn]*change;
@@ -1193 +1275 @@
         }
         if (upperBound(iColumn)) {
-          double change =upper_[iColumn]-oldPrimal-deltaX_[iColumn]-upperSlack_[iColumn];
+          CoinWorkDouble change =upper_[iColumn]-oldPrimal-deltaX_[iColumn]-upperSlack_[iColumn];
           c += upperSlack_[iColumn]*wVec_[iColumn];
           b += upperSlack_[iColumn]*deltaW_[iColumn]+wVec_[iColumn]*change;
@@ -1211 +1293 @@
     // maybe use inf and do line search
     // To check see if matches at current step
-    double step=CoinMin(actualPrimalStep_,actualDualStep_);
-    double next = c+b*step+a*step*step;
+    CoinWorkDouble step=CoinMin(actualPrimalStep_,actualDualStep_);
+    CoinWorkDouble next = c+b*step+a*step*step;
 #ifdef SOME_DEBUG
     printf("lin %g %g %g -> %g\n",a,b,c,
@@ -1230 +1312 @@
     a=0.0;
     b=0.0;
-    double ratio = actualDualStep_/actualPrimalStep_;
+    CoinWorkDouble ratio = actualDualStep_/actualPrimalStep_;
     for (iColumn=0;iColumn<numberTotal;iColumn++) {
-      double oldPrimal = solution_[iColumn];
+      CoinWorkDouble oldPrimal = solution_[iColumn];
       if (!flagged(iColumn)) {
         if (lowerBound(iColumn)) {
-          double change =oldPrimal+deltaX_[iColumn]-lowerSlack_[iColumn]-lower_[iColumn];
+          CoinWorkDouble change =oldPrimal+deltaX_[iColumn]-lowerSlack_[iColumn]-lower_[iColumn];
           b += lowerSlack_[iColumn]*deltaZ_[iColumn]*ratio+zVec_[iColumn]*change;
           a += deltaZ_[iColumn]*change*ratio;
         }
         if (upperBound(iColumn)) {
-          double change =upper_[iColumn]-oldPrimal-deltaX_[iColumn]-upperSlack_[iColumn];
+          CoinWorkDouble change =upper_[iColumn]-oldPrimal-deltaX_[iColumn]-upperSlack_[iColumn];
           b += upperSlack_[iColumn]*deltaW_[iColumn]*ratio+wVec_[iColumn]*change;
           a += deltaW_[iColumn]*change*ratio;
@@ -1250 +1332 @@
     // To check see if matches at current step
     step=actualPrimalStep_;
-    double next2 = c+b*step+a*step*step;
+    CoinWorkDouble next2 = c+b*step+a*step*step;
     if (next2>next) {
       actualPrimalStep_=CoinMin(actualPrimalStep_,actualDualStep_);
@@ -1275 +1357 @@
 #ifdef FULL_DEBUG
   if (phase==3){
-    double minBeta = 0.1*mu_;
-    double maxBeta = 10.0*mu_;
+    CoinWorkDouble minBeta = 0.1*mu_;
+    CoinWorkDouble maxBeta = 10.0*mu_;
     for (iColumn=0;iColumn<numberRows_+numberColumns_;iColumn++) {
       if (!flagged(iColumn)) {
         if (lowerBound(iColumn)) {
-          double change = -rhsL_[iColumn] + deltaX_[iColumn];
-          double dualValue=zVec_[iColumn]+actualDualStep_*deltaZ_[iColumn];
-          double primalValue=lowerSlack_[iColumn]+actualPrimalStep_*change;
-          double gapProduct=dualValue*primalValue;
+          CoinWorkDouble change = -rhsL_[iColumn] + deltaX_[iColumn];
+          CoinWorkDouble dualValue=zVec_[iColumn]+actualDualStep_*deltaZ_[iColumn];
+          CoinWorkDouble primalValue=lowerSlack_[iColumn]+actualPrimalStep_*change;
+          CoinWorkDouble gapProduct=dualValue*primalValue;
           if (delta2Z_[iColumn]<minBeta||delta2Z_[iColumn]>maxBeta)
             printf("3lower %d primal %g, dual %g, gap %g, old gap %g\n",
@@ -1289 +1371 @@
         }  
         if (upperBound(iColumn)) {
-          double change = rhsU_[iColumn]-deltaX_[iColumn];
-          double dualValue=wVec_[iColumn]+actualDualStep_*deltaW_[iColumn];
-          double primalValue=upperSlack_[iColumn]+actualPrimalStep_*change;
-          double gapProduct=dualValue*primalValue;
+          CoinWorkDouble change = rhsU_[iColumn]-deltaX_[iColumn];
+          CoinWorkDouble dualValue=wVec_[iColumn]+actualDualStep_*deltaW_[iColumn];
+          CoinWorkDouble primalValue=upperSlack_[iColumn]+actualPrimalStep_*change;
+          CoinWorkDouble gapProduct=dualValue*primalValue;
           if (delta2W_[iColumn]<minBeta||delta2W_[iColumn]>maxBeta)
             printf("3upper %d primal %g, dual %g, gap %g, old gap %g\n",
@@ -1303 +1385 @@
 #ifdef SOME_DEBUG_not
   {
-    double largestL=0.0;
-    double smallestL=COIN_DBL_MAX;
-    double largestU=0.0;
-    double smallestU=COIN_DBL_MAX;
-    double sumL=0.0;
-    double sumU=0.0;
+    CoinWorkDouble largestL=0.0;
+    CoinWorkDouble smallestL=COIN_DBL_MAX;
+    CoinWorkDouble largestU=0.0;
+    CoinWorkDouble smallestU=COIN_DBL_MAX;
+    CoinWorkDouble sumL=0.0;
+    CoinWorkDouble sumU=0.0;
     int nL=0;
     int nU=0;
@@ -1314 +1396 @@
       if (!flagged(iColumn)) {
         if (lowerBound(iColumn)) {
-          double change = -rhsL_[iColumn] + deltaX_[iColumn];
-          double dualValue=zVec_[iColumn]+actualDualStep_*deltaZ_[iColumn];
-          double primalValue=lowerSlack_[iColumn]+actualPrimalStep_*change;
-          double gapProduct=dualValue*primalValue;
+          CoinWorkDouble change = -rhsL_[iColumn] + deltaX_[iColumn];
+          CoinWorkDouble dualValue=zVec_[iColumn]+actualDualStep_*deltaZ_[iColumn];
+          CoinWorkDouble primalValue=lowerSlack_[iColumn]+actualPrimalStep_*change;
+          CoinWorkDouble gapProduct=dualValue*primalValue;
           largestL = CoinMax(largestL,gapProduct);
           smallestL = CoinMin(smallestL,gapProduct);
@@ -1324 +1406 @@
         }  
         if (upperBound(iColumn)) {
-          double change = rhsU_[iColumn]-deltaX_[iColumn];
-          double dualValue=wVec_[iColumn]+actualDualStep_*deltaW_[iColumn];
-          double primalValue=upperSlack_[iColumn]+actualPrimalStep_*change;
-          double gapProduct=dualValue*primalValue;
+          CoinWorkDouble change = rhsU_[iColumn]-deltaX_[iColumn];
+          CoinWorkDouble dualValue=wVec_[iColumn]+actualDualStep_*deltaW_[iColumn];
+          CoinWorkDouble primalValue=upperSlack_[iColumn]+actualPrimalStep_*change;
+          CoinWorkDouble gapProduct=dualValue*primalValue;
           largestU = CoinMax(largestU,gapProduct);
           smallestU = CoinMin(smallestU,gapProduct);
@@ -1335 +1417 @@
       } 
     }
-    double mu = (sumL+sumU)/(static_cast<double> (nL+nU));
+    CoinWorkDouble mu = (sumL+sumU)/(static_cast<CoinWorkDouble> (nL+nU));
 
-    double minBeta = 0.1*mu;
-    double maxBeta = 10.0*mu;
+    CoinWorkDouble minBeta = 0.1*mu;
+    CoinWorkDouble maxBeta = 10.0*mu;
     int nBL=0;
     int nAL=0;
@@ -1346 +1428 @@
       if (!flagged(iColumn)) {
         if (lowerBound(iColumn)) {
-          double change = -rhsL_[iColumn] + deltaX_[iColumn];
-          double dualValue=zVec_[iColumn]+actualDualStep_*deltaZ_[iColumn];
-          double primalValue=lowerSlack_[iColumn]+actualPrimalStep_*change;
-          double gapProduct=dualValue*primalValue;
+          CoinWorkDouble change = -rhsL_[iColumn] + deltaX_[iColumn];
+          CoinWorkDouble dualValue=zVec_[iColumn]+actualDualStep_*deltaZ_[iColumn];
+          CoinWorkDouble primalValue=lowerSlack_[iColumn]+actualPrimalStep_*change;
+          CoinWorkDouble gapProduct=dualValue*primalValue;
           if (gapProduct<minBeta)
             nBL++;
@@ -1359 +1441 @@
         }  
         if (upperBound(iColumn)) {
-          double change = rhsU_[iColumn]-deltaX_[iColumn];
-          double dualValue=wVec_[iColumn]+actualDualStep_*deltaW_[iColumn];
-          double primalValue=upperSlack_[iColumn]+actualPrimalStep_*change;
-          double gapProduct=dualValue*primalValue;
+          CoinWorkDouble change = rhsU_[iColumn]-deltaX_[iColumn];
+          CoinWorkDouble dualValue=wVec_[iColumn]+actualDualStep_*deltaW_[iColumn];
+          CoinWorkDouble primalValue=upperSlack_[iColumn]+actualPrimalStep_*change;
+          CoinWorkDouble gapProduct=dualValue*primalValue;
           if (gapProduct<minBeta)
             nBU++;
@@ -1384 +1466 @@
 /* Does solve. region1 is for deltaX (columns+rows), region2 for deltaPi (rows) */
 void 
-ClpPredictorCorrector::solveSystem(double * region1, double * region2,
-                                   const double * region1In, const double * region2In,
-                                   const double * saveRegion1, const double * saveRegion2,
+ClpPredictorCorrector::solveSystem(CoinWorkDouble * region1, CoinWorkDouble * region2,
+                                   const CoinWorkDouble * region1In, const CoinWorkDouble * region2In,
+                                   const CoinWorkDouble * saveRegion1, const CoinWorkDouble * saveRegion2,
                                    bool gentleRefine)
 {
@@ -1406 +1488 @@
     multiplyAdd(region1+numberColumns_,numberRows_,-1.0,region2,1.0);
     matrix_->times(1.0,region1,region2);
-    double maximumRHS = maximumAbsElement(region2,numberRows_);
-    double scale=1.0;
-    double unscale=1.0;
+    CoinWorkDouble maximumRHS = maximumAbsElement(region2,numberRows_);
+    CoinWorkDouble scale=1.0;
+    CoinWorkDouble unscale=1.0;
     if (maximumRHS>1.0e-30) {
       if (maximumRHS<=0.5) {
-        double factor=2.0;
+        CoinWorkDouble factor=2.0;
         while (maximumRHS<=0.5) {
           maximumRHS*=factor;
@@ -1417 +1499 @@
         } /* endwhile */
       } else if (maximumRHS>=2.0&&maximumRHS<=COIN_DBL_MAX) {
-        double factor=0.5;
+        CoinWorkDouble factor=0.5;
         while (maximumRHS>=2.0) {
           maximumRHS*=factor;
@@ -1444 +1526 @@
   if (saveRegion2) {
     //refine?
-    double scaleX=1.0;
+    CoinWorkDouble scaleX=1.0;
     if (gentleRefine) 
       scaleX=0.8;
@@ -1453 +1535 @@
 }
 // findDirectionVector.
-double ClpPredictorCorrector::findDirectionVector(const int phase)
+CoinWorkDouble ClpPredictorCorrector::findDirectionVector(const int phase)
 {
-  double projectionTolerance=projectionTolerance_;
+  CoinWorkDouble projectionTolerance=projectionTolerance_;
   //temporary
   //projectionTolerance=1.0e-15;
-  double errorCheck=0.9*maximumRHSError_/solutionNorm_;
+  CoinWorkDouble errorCheck=0.9*maximumRHSError_/solutionNorm_;
   if (errorCheck>primalTolerance()) {
     if (errorCheck<projectionTolerance) {
@@ -1468 +1550 @@
     } 
   } 
-  double * newError = new double [numberRows_];
+  CoinWorkDouble * newError = new CoinWorkDouble [numberRows_];
   int numberTotal = numberRows_+numberColumns_;
   //if flagged then entries zero so can do
   // For KKT separate out
-  double * region1Save=NULL;//for refinement
+  CoinWorkDouble * region1Save=NULL;//for refinement
   int iColumn;
   if (cholesky_->type()<20) {
@@ -1493 +1575 @@
   }
   bool goodSolve=false;
-  double * regionSave=NULL;//for refinement
+  CoinWorkDouble * regionSave=NULL;//for refinement
   int numberTries=0;
-  double relativeError=COIN_DBL_MAX;
-  double tryError=1.0e31;
+  CoinWorkDouble relativeError=COIN_DBL_MAX;
+  CoinWorkDouble tryError=1.0e31;
+  CoinWorkDouble saveMaximum=0.0;
+  double firstError=0.0;
+  double lastError=0.0;
   while (!goodSolve&&numberTries<30) {
-    double lastError=relativeError;
+    CoinWorkDouble lastError=relativeError;
     goodSolve=true;
-    double maximumRHS;
-    double saveMaximum;
+    CoinWorkDouble maximumRHS;
     maximumRHS = CoinMax(maximumAbsElement(deltaY_,numberRows_),1.0e-12);
-    saveMaximum = maximumRHS;
+    if (!numberTries)
+      saveMaximum = maximumRHS;
     if (cholesky_->type()<20) {
       // no kkt
-      double scale=1.0;
-      double unscale=1.0;
+      CoinWorkDouble scale=1.0;
+      CoinWorkDouble unscale=1.0;
       if (maximumRHS>1.0e-30) {
         if (maximumRHS<=0.5) {
-          double factor=2.0;
+          CoinWorkDouble factor=2.0;
           while (maximumRHS<=0.5) {
             maximumRHS*=factor;
@@ -1516 +1601 @@
           } /* endwhile */
         } else if (maximumRHS>=2.0&&maximumRHS<=COIN_DBL_MAX) {
-          double factor=0.5;
+          CoinWorkDouble factor=0.5;
           while (maximumRHS>=2.0) {
             maximumRHS*=factor;
@@ -1544 +1629 @@
       if (numberTries) {
         //refine?
-        double scaleX=1.0;
+        CoinWorkDouble scaleX=1.0;
         if (lastError>1.0e-5) 
           scaleX=0.8;
@@ -1567 +1652 @@
     
     //now add in old Ax - doing extra checking
-    double maximumRHSError=0.0;
-    double maximumRHSChange=0.0;
+    CoinWorkDouble maximumRHSError=0.0;
+    CoinWorkDouble maximumRHSChange=0.0;
     int iRow;
     char * dropped = cholesky_->rowsDropped();
     for (iRow=0;iRow<numberRows_;iRow++) {
       if (!dropped[iRow]) {
-        double newValue=newError[iRow];
-        double oldValue=errorRegion_[iRow];
+        CoinWorkDouble newValue=newError[iRow];
+        CoinWorkDouble oldValue=errorRegion_[iRow];
         //severity of errors depend on signs
         //**later                                                             */
-        if (fabs(newValue)>maximumRHSChange) {
-          maximumRHSChange=fabs(newValue);
+        if (CoinAbs(newValue)>maximumRHSChange) {
+          maximumRHSChange=CoinAbs(newValue);
         } 
-        double result=newValue+oldValue;
-        if (fabs(result)>maximumRHSError) {
-          maximumRHSError=fabs(result);
+        CoinWorkDouble result=newValue+oldValue;
+        if (CoinAbs(result)>maximumRHSError) {
+          maximumRHSError=CoinAbs(result);
         } 
         newError[iRow]=result;
       } else {
-        double newValue=newError[iRow];
-        double oldValue=errorRegion_[iRow];
-        if (fabs(newValue)>maximumRHSChange) {
-          maximumRHSChange=fabs(newValue);
+        CoinWorkDouble newValue=newError[iRow];
+        CoinWorkDouble oldValue=errorRegion_[iRow];
+        if (CoinAbs(newValue)>maximumRHSChange) {
+          maximumRHSChange=CoinAbs(newValue);
         } 
-        double result=newValue+oldValue;
+        CoinWorkDouble result=newValue+oldValue;
         newError[iRow]=result;
         //newError[iRow]=0.0;
@@ -1602 +1687 @@
     if (relativeError>tryError) 
       relativeError=tryError;
+    if (numberTries==1)
+      firstError=relativeError;
     if (relativeError<lastError) {
+      lastError=relativeError;
       maximumRHSChange_= maximumRHSChange;
-      if (relativeError>1.0e-9
-          ||numberTries>1) {
-        handler_->message(CLP_BARRIER_ACCURACY,messages_)
-          <<phase<<numberTries<<relativeError
-          <<CoinMessageEol;
-      } 
       if (relativeError>projectionTolerance&&numberTries<=3) {
         //try and refine
@@ -1617 +1699 @@
       if (!goodSolve) {
         if (!regionSave) {
-          regionSave = new double [numberRows_];
+          regionSave = new CoinWorkDouble [numberRows_];
           if (cholesky_->type()>=20) 
-            region1Save = new double [numberTotal];
+            region1Save = new CoinWorkDouble [numberTotal];
         } 
         CoinMemcpyN(deltaY_,numberRows_,regionSave);
@@ -1653 +1735 @@
       } else {
         // disaster
-        CoinFillN(deltaX_,numberTotal,1.0);
-        CoinFillN(deltaY_,numberRows_,1.0);
+        CoinFillN(deltaX_,numberTotal,static_cast<CoinWorkDouble>(1.0));
+        CoinFillN(deltaY_,numberRows_,static_cast<CoinWorkDouble>(1.0));
         printf("bad cholesky\n");
       }
     }
   } /* endwhile */
+  if (firstError>1.0e-9||numberTries>1) {
+    handler_->message(CLP_BARRIER_ACCURACY,messages_)
+      <<phase<<numberTries<<static_cast<double>(firstError)
+      <<static_cast<double>(lastError)
+      <<CoinMessageEol;
+  } 
   delete [] regionSave;
   delete [] region1Save;
   delete [] newError;
   // now rest
-  double extra=eExtra;
+  CoinWorkDouble extra=eExtra;
   //multiplyAdd(deltaY_,numberRows_,1.0,deltaW_+numberColumns_,0.0);
   //CoinZeroN(deltaW_,numberColumns_);
@@ -1672 +1760 @@
     deltaSL_[iColumn]=0.0;
     deltaZ_[iColumn]=0.0;
-    double dd=deltaW_[iColumn];
+    CoinWorkDouble dd=deltaW_[iColumn];
     deltaW_[iColumn]=0.0;
     if (!flagged(iColumn)) {
-      double deltaX = deltaX_[iColumn];
+      CoinWorkDouble deltaX = deltaX_[iColumn];
       if (lowerBound(iColumn)) {
-        double zValue = rhsZ_[iColumn];
-        double gHat = zValue + zVec_[iColumn]*rhsL_[iColumn];
-        double slack = lowerSlack_[iColumn]+extra;
+        CoinWorkDouble zValue = rhsZ_[iColumn];
+        CoinWorkDouble gHat = zValue + zVec_[iColumn]*rhsL_[iColumn];
+        CoinWorkDouble slack = lowerSlack_[iColumn]+extra;
         deltaSL_[iColumn] = -rhsL_[iColumn]+deltaX;
         deltaZ_[iColumn]=(gHat-zVec_[iColumn]*deltaX)/slack;
       } 
       if (upperBound(iColumn)) {
-        double wValue = rhsW_[iColumn];
-        double hHat = wValue - wVec_[iColumn]*rhsU_[iColumn];
-        double slack = upperSlack_[iColumn]+extra;
+        CoinWorkDouble wValue = rhsW_[iColumn];
+        CoinWorkDouble hHat = wValue - wVec_[iColumn]*rhsU_[iColumn];
+        CoinWorkDouble slack = upperSlack_[iColumn]+extra;
         deltaSU_[iColumn] = rhsU_[iColumn]-deltaX;
         deltaW_[iColumn]=(hHat+wVec_[iColumn]*deltaX)/slack;
@@ -1692 +1780 @@
       if (0) {
         // different way of calculating
-        double gamma2 = gamma_*gamma_;
-        double dZ=0.0;
-        double dW=0.0;
-        double zValue = rhsZ_[iColumn];
-        double gHat = zValue + zVec_[iColumn]*rhsL_[iColumn];
-        double slackL = lowerSlack_[iColumn]+extra;
-        double wValue = rhsW_[iColumn];
-        double hHat = wValue - wVec_[iColumn]*rhsU_[iColumn];
-        double slackU = upperSlack_[iColumn]+extra;
-        double q = rhsC_[iColumn]+gamma2 * deltaX +dd;
+        CoinWorkDouble gamma2 = gamma_*gamma_;
+        CoinWorkDouble dZ=0.0;
+        CoinWorkDouble dW=0.0;
+        CoinWorkDouble zValue = rhsZ_[iColumn];
+        CoinWorkDouble gHat = zValue + zVec_[iColumn]*rhsL_[iColumn];
+        CoinWorkDouble slackL = lowerSlack_[iColumn]+extra;
+        CoinWorkDouble wValue = rhsW_[iColumn];
+        CoinWorkDouble hHat = wValue - wVec_[iColumn]*rhsU_[iColumn];
+        CoinWorkDouble slackU = upperSlack_[iColumn]+extra;
+        CoinWorkDouble q = rhsC_[iColumn]+gamma2 * deltaX +dd;
         if (primalR_)
           q += deltaX*primalR_[iColumn];
@@ -1728 +1816 @@
   }
 #if 0
-  double * check = new double[numberTotal];  
+  CoinWorkDouble * check = new CoinWorkDouble[numberTotal];  
   // Check out rhsC_
   multiplyAdd(deltaY_,numberRows_,-1.0,check+numberColumns_,0.0);
@@ -1736 +1824 @@
   for (iColumn=0;iColumn<numberTotal;iColumn++) {
     check[iColumn] += deltaZ_[iColumn]-deltaW_[iColumn];
-    if (fabs(check[iColumn]-rhsC_[iColumn])>1.0e-3)
+    if (CoinAbs(check[iColumn]-rhsC_[iColumn])>1.0e-3)
       printf("rhsC %d %g %g\n",iColumn,check[iColumn],rhsC_[iColumn]);
   }
@@ -1743 +1831 @@
     check[iColumn] += lowerSlack_[iColumn]*deltaZ_[iColumn]+
       zVec_[iColumn]*deltaSL_[iColumn];
-    if (fabs(check[iColumn]-rhsZ_[iColumn])>1.0e-3)
+    if (CoinAbs(check[iColumn]-rhsZ_[iColumn])>1.0e-3)
       printf("rhsZ %d %g %g\n",iColumn,check[iColumn],rhsZ_[iColumn]);
   }
@@ -1750 +1838 @@
     check[iColumn] += upperSlack_[iColumn]*deltaW_[iColumn]+
       wVec_[iColumn]*deltaSU_[iColumn];
-    if (fabs(check[iColumn]-rhsW_[iColumn])>1.0e-3)
+    if (CoinAbs(check[iColumn]-rhsW_[iColumn])>1.0e-3)
       printf("rhsW %d %g %g\n",iColumn,check[iColumn],rhsW_[iColumn]);
   }
@@ -1762 +1850 @@
   int numberTotal = numberRows_+numberColumns_;
   int iColumn;
-  double tolerance = primalTolerance();
+  CoinWorkDouble tolerance = primalTolerance();
   // See if quadratic objective
 #ifndef NO_RTTI
@@ -1783 +1871 @@
       clearFixed(iColumn);
   }
-  double initialValue =1.0e-12;
 
-  double maximumObjective=0.0;
-  double objectiveNorm2=0.0;
+  CoinWorkDouble maximumObjective=0.0;
+  CoinWorkDouble objectiveNorm2=0.0;
   getNorms(cost_,numberTotal,maximumObjective,objectiveNorm2);
   if (!maximumObjective) {
     maximumObjective=1.0; // objective all zero
   } 
-  objectiveNorm2=sqrt(objectiveNorm2)/static_cast<double> (numberTotal);
+  objectiveNorm2=CoinSqrt(objectiveNorm2)/static_cast<CoinWorkDouble> (numberTotal);
   objectiveNorm_=maximumObjective;
   scaleFactor_=1.0;
@@ -1809 +1896 @@
   if (quadraticObj) {
     // If scaled then really scale matrix
-    double scaleFactor = 
+    CoinWorkDouble scaleFactor = 
       scaleFactor_*optimizationDirection_*objectiveScale_*
       rhsScale_;
@@ -1819 +1906 @@
       double * quadraticElement = quadratic->getMutableElements();
       int numberColumns = quadratic->getNumCols();
-      double scale = 1.0/scaleFactor;
+      CoinWorkDouble scale = 1.0/scaleFactor;
       if (scalingFlag_>0&&rowScale_) {
         for (int iColumn=0;iColumn<numberColumns;iColumn++) {
-          double scaleI = columnScale_[iColumn]*scale;
+          CoinWorkDouble scaleI = columnScale_[iColumn]*scale;
           for (CoinBigIndex j=columnQuadraticStart[iColumn];
                j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
             int jColumn = columnQuadratic[j];
-            double scaleJ = columnScale_[jColumn];
+            CoinWorkDouble scaleJ = columnScale_[jColumn];
             quadraticElement[j] *= scaleI*scaleJ;
-            objectiveNorm_ = CoinMax(objectiveNorm_,fabs(quadraticElement[j]));
+            objectiveNorm_ = CoinMax(objectiveNorm_,CoinAbs(quadraticElement[j]));
           }
         }
@@ -1837 +1924 @@
                j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
             quadraticElement[j] *= scale;
-            objectiveNorm_ = CoinMax(objectiveNorm_,fabs(quadraticElement[j]));
+            objectiveNorm_ = CoinMax(objectiveNorm_,CoinAbs(quadraticElement[j]));
           }
         }
@@ -1848 +1935 @@
   //DZ in lowerDual
   //DW in upperDual
-  double infiniteCheck=1.0e40;
-  //double     fakeCheck=1.0e10;
+  CoinWorkDouble infiniteCheck=1.0e40;
+  //CoinWorkDouble     fakeCheck=1.0e10;
   //use deltaX region for work region
   for (iColumn=0;iColumn<numberTotal;iColumn++) {
-    double primalValue = solution_[iColumn];
+    CoinWorkDouble primalValue = solution_[iColumn];
     clearFlagged(iColumn);
     clearFixedOrFree(iColumn);
@@ -1863 +1950 @@
       diagonal_[iColumn]=1.0;
       deltaX_[iColumn]=1.0;
-      double lowerValue=lower_[iColumn];
-      double upperValue=upper_[iColumn];
+      CoinWorkDouble lowerValue=lower_[iColumn];
+      CoinWorkDouble upperValue=upper_[iColumn];
       if (lowerValue>-infiniteCheck) {
         if (upperValue<infiniteCheck) {
@@ -1953 +2040 @@
     solveSystem(deltaX_,errorRegion_,solution_,NULL,NULL,NULL,false);
   }
-  initialValue=1.0e2;
+  CoinWorkDouble initialValue = 1.0e2;
   if (rhsNorm_*1.0e-2>initialValue) {
     initialValue=rhsNorm_*1.0e-2;
   } 
   //initialValue = CoinMax(1.0,rhsNorm_);
-  double smallestBoundDifference=COIN_DBL_MAX;
-  double * fakeSolution = deltaX_;
+  CoinWorkDouble smallestBoundDifference=COIN_DBL_MAX;
+  CoinWorkDouble * fakeSolution = deltaX_;
   for ( iColumn=0;iColumn<numberTotal;iColumn++) {
     if (!flagged(iColumn)) {
@@ -1975 +2062 @@
   solutionNorm_=1.0e-12;
   handler_->message(CLP_BARRIER_SAFE,messages_)
-    <<initialValue<<objectiveNorm_
+    <<static_cast<double>(initialValue)<<static_cast<double>(objectiveNorm_)
     <<CoinMessageEol;
-  double extra=1.0e-10;
-  double largeGap=1.0e15;
-  //double safeObjectiveValue=2.0*objectiveNorm_;
-  double safeObjectiveValue=objectiveNorm_+1.0;
-  double safeFree=1.0e-1*initialValue;
+  CoinWorkDouble extra=1.0e-10;
+  CoinWorkDouble largeGap=1.0e15;
+  //CoinWorkDouble safeObjectiveValue=2.0*objectiveNorm_;
+  CoinWorkDouble safeObjectiveValue=objectiveNorm_+1.0;
+  CoinWorkDouble safeFree=1.0e-1*initialValue;
   //printf("normal safe dual value of %g, primal value of %g\n",
   // safeObjectiveValue,initialValue);
@@ -1988 +2075 @@
   //printf("temp safe dual value of %g, primal value of %g\n",
   // safeObjectiveValue,initialValue);
-  double zwLarge=1.0e2*initialValue;
+  CoinWorkDouble zwLarge=1.0e2*initialValue;
   //zwLarge=1.0e40;
   if (cholesky_->choleskyCondition()<0.0&&cholesky_->type()<20) {
@@ -1996 +2083 @@
     safeFree *= 10.0;
   }
-  double gamma2 = gamma_*gamma_; // gamma*gamma will be added to diagonal
-#if 0
-  fakeSolution[0 ] =   0.072310129 ;
-  fakeSolution[1 ] =   0.053083871; 
-  fakeSolution[2 ] =      0.178127; 
-  fakeSolution[3 ] =    0.13215151; 
-  fakeSolution[4 ] =   0.072715642; 
-  fakeSolution[5 ] =    0.15680727; 
-  fakeSolution[6 ] =    0.16841689; 
-  fakeSolution[7 ] =   0.093612798 ;
-  fakeSolution[8 ] =   0.072774891 ;
-  fakeSolution[9]=1.0;
-  initialValue=1.0e-5;
-  safeObjectiveValue=1.0e-5;
-#endif
+  CoinWorkDouble gamma2 = gamma_*gamma_; // gamma*gamma will be added to diagonal
   // First do primal side
   for ( iColumn=0;iColumn<numberTotal;iColumn++) {
     if (!flagged(iColumn)) {
-      double lowerValue=lower_[iColumn];
-      double upperValue=upper_[iColumn];
-      double newValue;
-      double setPrimal=initialValue;
+      CoinWorkDouble lowerValue=lower_[iColumn];
+      CoinWorkDouble upperValue=upper_[iColumn];
+      CoinWorkDouble newValue;
+      CoinWorkDouble setPrimal=initialValue;
       if (quadraticObj) {
         // perturb primal solution a bit
@@ -2026 +2099 @@
           //upper and lower bounds
           if (upperValue-lowerValue>2.0*setPrimal) {
-            double fakeValue=fakeSolution[iColumn];
+            CoinWorkDouble fakeValue=fakeSolution[iColumn];
             if (fakeValue<lowerValue+setPrimal) {
               fakeValue=lowerValue+setPrimal;
@@ -2039 +2112 @@
         } else {
           //just lower bound
-          double fakeValue=fakeSolution[iColumn];
+          CoinWorkDouble fakeValue=fakeSolution[iColumn];
           if (fakeValue<lowerValue+setPrimal) {
             fakeValue=lowerValue+setPrimal;
@@ -2048 +2121 @@
         if (upperBound(iColumn)) {
           //just upper bound
-          double fakeValue=fakeSolution[iColumn];
+          CoinWorkDouble fakeValue=fakeSolution[iColumn];
           if (fakeValue>upperValue-setPrimal) {
             fakeValue=upperValue-setPrimal;
@@ -2093 +2166 @@
     quadraticElement = quadratic->getElements();
   }
-  double quadraticNorm=0.0;
+  CoinWorkDouble quadraticNorm=0.0;
   for ( iColumn=0;iColumn<numberTotal;iColumn++) {
     if (!flagged(iColumn)) {
-      double primalValue=solution_[iColumn];
-      double lowerValue=lower_[iColumn];
-      double upperValue=upper_[iColumn];
+      CoinWorkDouble primalValue=solution_[iColumn];
+      CoinWorkDouble lowerValue=lower_[iColumn];
+      CoinWorkDouble upperValue=upper_[iColumn];
       // Do dj
-      double reducedCost=cost_[iColumn];
+      CoinWorkDouble reducedCost=cost_[iColumn];
       if (lowerBound(iColumn)) {
         reducedCost+=linearPerturbation_;
@@ -2111 +2184 @@
              j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
           int jColumn = columnQuadratic[j];
-          double valueJ = solution_[jColumn];
-          double elementValue = quadraticElement[j];
+          CoinWorkDouble valueJ = solution_[jColumn];
+          CoinWorkDouble elementValue = quadraticElement[j];
           reducedCost += valueJ*elementValue;
         }
-        quadraticNorm = CoinMax(quadraticNorm,fabs(reducedCost));
+        quadraticNorm = CoinMax(quadraticNorm,CoinAbs(reducedCost));
       }
       dj_[iColumn]=reducedCost;
@@ -2132 +2205 @@
   for ( iColumn=0;iColumn<numberTotal;iColumn++) {
     if (!flagged(iColumn)) {
-      double primalValue=solution_[iColumn];
-      double lowerValue=lower_[iColumn];
-      double upperValue=upper_[iColumn];
-      double reducedCost=dj_[iColumn];
-      double low=0.0;
-      double high=0.0;
+      CoinWorkDouble primalValue=solution_[iColumn];
+      CoinWorkDouble lowerValue=lower_[iColumn];
+      CoinWorkDouble upperValue=upper_[iColumn];
+      CoinWorkDouble reducedCost=dj_[iColumn];
+      CoinWorkDouble low=0.0;
+      CoinWorkDouble high=0.0;
       if (lowerBound(iColumn)) {
         if (upperBound(iColumn)) {
@@ -2148 +2221 @@
             high=initialValue;
           }
-          double s = low+extra;
-          double ratioZ;
+          CoinWorkDouble s = low+extra;
+          CoinWorkDouble ratioZ;
           if (s<zwLarge) {
             ratioZ=1.0;
           } else {
-            ratioZ=sqrt(zwLarge/s);
+            ratioZ=CoinSqrt(zwLarge/s);
           } 
-          double t = high+extra;
-          double ratioT;
+          CoinWorkDouble t = high+extra;
+          CoinWorkDouble ratioT;
           if (t<zwLarge) {
             ratioT=1.0;
           } else {
-            ratioT=sqrt(zwLarge/t);
+            ratioT=CoinSqrt(zwLarge/t);
           } 
           //modify s and t
@@ -2179 +2252 @@
             wVec_[iColumn]=CoinMax(-reducedCost , safeObjectiveValue*ratioT);
           }
-          double gammaTerm = gamma2;
+          CoinWorkDouble gammaTerm = gamma2;
           if (primalR_)
             gammaTerm += primalR_[iColumn];
@@ -2188 +2261 @@
           low=primalValue-lowerValue;
           high=0.0;
-          double s = low+extra;
-          double ratioZ;
+          CoinWorkDouble s = low+extra;
+          CoinWorkDouble ratioZ;
           if (s<zwLarge) {
             ratioZ=1.0;
           } else {
-            ratioZ=sqrt(zwLarge/s);
+            ratioZ=CoinSqrt(zwLarge/s);
           } 
           //modify s
@@ -2207 +2280 @@
             wVec_[iColumn]=0.0;
           } 
-          double gammaTerm = gamma2;
+          CoinWorkDouble gammaTerm = gamma2;
           if (primalR_)
             gammaTerm += primalR_[iColumn];
@@ -2217 +2290 @@
           low=0.0;
           high=upperValue-primalValue;
-          double t = high+extra;
-          double ratioT;
+          CoinWorkDouble t = high+extra;
+          CoinWorkDouble ratioT;
           if (t<zwLarge) {
             ratioT=1.0;
           } else {
-            ratioT=sqrt(zwLarge/t);
+            ratioT=CoinSqrt(zwLarge/t);
           } 
           //modify t
@@ -2236 +2309 @@
             wVec_[iColumn]=CoinMax(-reducedCost , safeObjectiveValue*ratioT);
           } 
-          double gammaTerm = gamma2;
+          CoinWorkDouble gammaTerm = gamma2;
           if (primalR_)
             gammaTerm += primalR_[iColumn];
@@ -2249 +2322 @@
   if (solution_[0]>0.0) {
     for (int i=0;i<numberTotal;i++)
-      printf("%d %.18g %.18g %.18g %.18g %.18g %.18g %.18g\n",i,fabs(solution_[i]),
-             diagonal_[i],fabs(dj_[i]),
+      printf("%d %.18g %.18g %.18g %.18g %.18g %.18g %.18g\n",i,CoinAbs(solution_[i]),
+             diagonal_[i],CoinAbs(dj_[i]),
              lowerSlack_[i],zVec_[i],
              upperSlack_[i],wVec_[i]);
   } else {
     for (int i=0;i<numberTotal;i++)
-      printf("%d %.18g %.18g %.18g %.18g %.18g %.18g %.18g\n",i,fabs(solution_[i]),
-             diagonal_[i],fabs(dj_[i]),
+      printf("%d %.18g %.18g %.18g %.18g %.18g %.18g %.18g\n",i,CoinAbs(solution_[i]),
+             diagonal_[i],CoinAbs(dj_[i]),
              upperSlack_[i],wVec_[i],
              lowerSlack_[i],zVec_[i] );
@@ -2266 +2339 @@
 // complementarityGap.  Computes gap
 //phase 0=as is , 1 = after predictor , 2 after corrector
-double ClpPredictorCorrector::complementarityGap(int & numberComplementarityPairs,
+CoinWorkDouble ClpPredictorCorrector::complementarityGap(int & numberComplementarityPairs,
                                                  int & numberComplementarityItems,
                                                  const int phase)
 {
-  double gap=0.0;
+  CoinWorkDouble gap=0.0;
   //seems to be same coding for phase = 1 or 2
   numberComplementarityPairs=0;
   numberComplementarityItems=0;
   int numberTotal = numberRows_+numberColumns_;
-  double toleranceGap=0.0;
-  double largestGap=0.0;
-  double smallestGap=COIN_DBL_MAX;
+  CoinWorkDouble toleranceGap=0.0;
+  CoinWorkDouble largestGap=0.0;
+  CoinWorkDouble smallestGap=COIN_DBL_MAX;
   //seems to be same coding for phase = 1 or 2
   int numberNegativeGaps=0;
-  double sumNegativeGap=0.0;
-  double largeGap=1.0e2*solutionNorm_;
+  CoinWorkDouble sumNegativeGap=0.0;
+  CoinWorkDouble largeGap=1.0e2*solutionNorm_;
   if (largeGap<1.0e10) {
     largeGap=1.0e10;
   } 
   largeGap=1.0e30;
-  double dualTolerance =  dblParam_[ClpDualTolerance];
-  double primalTolerance =  dblParam_[ClpPrimalTolerance];
+  CoinWorkDouble dualTolerance =  dblParam_[ClpDualTolerance];
+  CoinWorkDouble primalTolerance =  dblParam_[ClpPrimalTolerance];
   dualTolerance=dualTolerance/scaleFactor_;
   for (int iColumn=0;iColumn<numberTotal;iColumn++) {
@@ -2293 +2366 @@
       numberComplementarityPairs++;
       //can collapse as if no lower bound both zVec and deltaZ 0.0
-      double newZ=0.0;
-      double newW=0.0;
+      CoinWorkDouble newZ=0.0;
+      CoinWorkDouble newW=0.0;
       if (lowerBound(iColumn)) {
         numberComplementarityItems++;
-        double dualValue;
-        double primalValue;
+        CoinWorkDouble dualValue;
+        CoinWorkDouble primalValue;
         if (!phase) {
           dualValue=zVec_[iColumn];
           primalValue=lowerSlack_[iColumn];
         } else {
-          double change;
+          CoinWorkDouble change;
           change =solution_[iColumn]+deltaX_[iColumn]-lowerSlack_[iColumn]-lower_[iColumn];
           dualValue=zVec_[iColumn]+actualDualStep_*deltaZ_[iColumn];
@@ -2313 +2386 @@
           primalValue=largeGap;
         } 
-        double gapProduct=dualValue*primalValue;
+        CoinWorkDouble gapProduct=dualValue*primalValue;
         if (gapProduct<0.0) {
           //cout<<"negative gap component "<<iColumn<<" "<<dualValue<<" "<<
@@ -2334 +2407 @@
       if (upperBound(iColumn)) {
         numberComplementarityItems++;
-        double dualValue;
-        double primalValue;
+        CoinWorkDouble dualValue;
+        CoinWorkDouble primalValue;
         if (!phase) {
           dualValue=wVec_[iColumn];
           primalValue=upperSlack_[iColumn];
         } else {
-          double change;
+          CoinWorkDouble change;
           change =upper_[iColumn]-solution_[iColumn]-deltaX_[iColumn]-upperSlack_[iColumn];
           dualValue=wVec_[iColumn]+actualDualStep_*deltaW_[iColumn];
@@ -2350 +2423 @@
           primalValue=largeGap;
         } 
-        double gapProduct=dualValue*primalValue;
+        CoinWorkDouble gapProduct=dualValue*primalValue;
         if (gapProduct<0.0) {
           //cout<<"negative gap component "<<iColumn<<" "<<dualValue<<" "<<
@@ -2374 +2447 @@
   if (!phase&&numberNegativeGaps) {
       handler_->message(CLP_BARRIER_NEGATIVE_GAPS,messages_)
-    <<numberNegativeGaps<<sumNegativeGap
+        <<numberNegativeGaps<<static_cast<double>(sumNegativeGap)
     <<CoinMessageEol;
   } 
@@ -2393 +2466 @@
 void ClpPredictorCorrector::setupForSolve(const int phase)
 {
-  double extra =eExtra;
+  CoinWorkDouble extra =eExtra;
   int numberTotal = numberRows_ + numberColumns_;
   int iColumn;
@@ -2399 +2472 @@
   printf("phase %d in setupForSolve, mu %.18g\n",phase,mu_);
 #endif
-  double gamma2 = gamma_*gamma_; // gamma*gamma will be added to diagonal
+  CoinWorkDouble gamma2 = gamma_*gamma_; // gamma*gamma will be added to diagonal
+  CoinWorkDouble * dualArray = reinterpret_cast<CoinWorkDouble *>(dual_);
   switch (phase) {
   case 0:
@@ -2405 +2479 @@
     if (delta_||dualR_) {
       // add in regularization
-      double delta2 = delta_*delta_;
+      CoinWorkDouble delta2 = delta_*delta_;
       for (int iRow=0;iRow<numberRows_;iRow++) {
-        rhsB_[iRow] -= delta2*dual_[iRow];
+        rhsB_[iRow] -= delta2*dualArray[iRow];
         if (dualR_)
-          rhsB_[iRow] -= dualR_[iRow]*dual_[iRow];
+          rhsB_[iRow] -= dualR_[iRow]*dualArray[iRow];
       }
     }
@@ -2445 +2519 @@
 #if 0
     for (int i=0;i<3;i++) {
-      if (!fabs(rhsZ_[i]))
+      if (!CoinAbs(rhsZ_[i]))
         rhsZ_[i]=0.0;
-      if (!fabs(rhsW_[i]))
+      if (!CoinAbs(rhsW_[i]))
         rhsW_[i]=0.0;
-      if (!fabs(rhsU_[i]))
+      if (!CoinAbs(rhsU_[i]))
         rhsU_[i]=0.0;
-      if (!fabs(rhsL_[i]))
+      if (!CoinAbs(rhsL_[i]))
         rhsL_[i]=0.0;
     }
@@ -2499 +2573 @@
 #if 0
     for (int i=0;i<numberTotal;i++) {
-      if (!fabs(rhsZ_[i]))
+      if (!CoinAbs(rhsZ_[i]))
         rhsZ_[i]=0.0;
-      if (!fabs(rhsW_[i]))
+      if (!CoinAbs(rhsW_[i]))
         rhsW_[i]=0.0;
-      if (!fabs(rhsU_[i]))
+      if (!CoinAbs(rhsU_[i]))
         rhsU_[i]=0.0;
-      if (!fabs(rhsL_[i]))
+      if (!CoinAbs(rhsL_[i]))
         rhsL_[i]=0.0;
     }
     if (solution_[0]>0.0) {
       for (int i=0;i<numberTotal;i++)
-        printf("%d %.18g %.18g %.18g %.18g %.18g %.18g %.18g\n",i,fabs(solution_[i]),
-               diagonal_[i],fabs(dj_[i]),
+        printf("%d %.18g %.18g %.18g %.18g %.18g %.18g %.18g\n",i,CoinAbs(solution_[i]),
+               diagonal_[i],CoinAbs(dj_[i]),
                lowerSlack_[i],zVec_[i],
                upperSlack_[i],wVec_[i]);
       for (int i=0;i<numberTotal;i++)
-        printf("%d %.18g %.18g %.18g %.18g %.18g\n",i,fabs(rhsC_[i]),
+        printf("%d %.18g %.18g %.18g %.18g %.18g\n",i,CoinAbs(rhsC_[i]),
                rhsZ_[i],rhsL_[i],
                rhsW_[i],rhsU_[i]);
     } else {
       for (int i=0;i<numberTotal;i++)
-        printf("%d %.18g %.18g %.18g %.18g %.18g %.18g %.18g\n",i,fabs(solution_[i]),
-               diagonal_[i],fabs(dj_[i]),
+        printf("%d %.18g %.18g %.18g %.18g %.18g %.18g %.18g\n",i,CoinAbs(solution_[i]),
+               diagonal_[i],CoinAbs(dj_[i]),
                upperSlack_[i],wVec_[i],
                lowerSlack_[i],zVec_[i] );
       for (int i=0;i<numberTotal;i++)
-        printf("%d %.18g %.18g %.18g %.18g %.18g\n",i,fabs(rhsC_[i]),
+        printf("%d %.18g %.18g %.18g %.18g %.18g\n",i,CoinAbs(rhsC_[i]),
                rhsW_[i],rhsU_[i],
                rhsZ_[i],rhsL_[i]);
@@ -2549 +2623 @@
   case 3:
     {
-      double minBeta = 0.1*mu_;
-      double maxBeta = 10.0*mu_;
-      double dualStep = CoinMin(1.0,actualDualStep_+0.1);
-      double primalStep = CoinMin(1.0,actualPrimalStep_+0.1);
+      CoinWorkDouble minBeta = 0.1*mu_;
+      CoinWorkDouble maxBeta = 10.0*mu_;
+      CoinWorkDouble dualStep = CoinMin(1.0,actualDualStep_+0.1);
+      CoinWorkDouble primalStep = CoinMin(1.0,actualPrimalStep_+0.1);
 #ifdef SOME_DEBUG
       printf("good complementarity range %g to %g\n",minBeta,maxBeta);
@@ -2561 +2635 @@
         if (!flagged(iColumn)) {
           if (lowerBound(iColumn)) {
-            double change = -rhsL_[iColumn] + deltaX_[iColumn];
-            double dualValue=zVec_[iColumn]+dualStep*deltaZ_[iColumn];
-            double primalValue=lowerSlack_[iColumn]+primalStep*change;
-            double gapProduct=dualValue*primalValue;
+            CoinWorkDouble change = -rhsL_[iColumn] + deltaX_[iColumn];
+            CoinWorkDouble dualValue=zVec_[iColumn]+dualStep*deltaZ_[iColumn];
+            CoinWorkDouble primalValue=lowerSlack_[iColumn]+primalStep*change;
+            CoinWorkDouble gapProduct=dualValue*primalValue;
             if (gapProduct>0.0&&dualValue<0.0)
               gapProduct = - gapProduct;
@@ -2573 +2647 @@
                      iColumn,primalValue,dualValue,gapProduct);
 #endif
-            double value= 0.0;
+            CoinWorkDouble value= 0.0;
             if (gapProduct<minBeta) {
               value= 2.0*(minBeta-gapProduct);
@@ -2586 +2660 @@
           }  
           if (upperBound(iColumn)) {
-            double change = rhsU_[iColumn]-deltaX_[iColumn];
-            double dualValue=wVec_[iColumn]+dualStep*deltaW_[iColumn];
-            double primalValue=upperSlack_[iColumn]+primalStep*change;
-            double gapProduct=dualValue*primalValue;
+            CoinWorkDouble change = rhsU_[iColumn]-deltaX_[iColumn];
+            CoinWorkDouble dualValue=wVec_[iColumn]+dualStep*deltaW_[iColumn];
+            CoinWorkDouble primalValue=upperSlack_[iColumn]+primalStep*change;
+            CoinWorkDouble gapProduct=dualValue*primalValue;
             if (gapProduct>0.0&&dualValue<0.0)
               gapProduct = - gapProduct;
@@ -2598 +2672 @@
                      iColumn,primalValue,dualValue,gapProduct);
 #endif
-            double value= 0.0;
+            CoinWorkDouble value= 0.0;
             if (gapProduct<minBeta) {
               value= (minBeta-gapProduct);
@@ -2615 +2689 @@
   if (cholesky_->type()<20) {
     for (iColumn=0;iColumn<numberTotal;iColumn++) {
-      double value = rhsC_[iColumn];
-      double zValue = rhsZ_[iColumn];
-      double wValue = rhsW_[iColumn];
+      CoinWorkDouble value = rhsC_[iColumn];
+      CoinWorkDouble zValue = rhsZ_[iColumn];
+      CoinWorkDouble wValue = rhsW_[iColumn];
 #if 0
 #if 1
@@ -2644 +2718 @@
 #else
       if (lowerBound(iColumn)) {
-        double gHat = zValue + zVec_[iColumn]*rhsL_[iColumn];
+        CoinWorkDouble gHat = zValue + zVec_[iColumn]*rhsL_[iColumn];
         value -= gHat/(lowerSlack_[iColumn]+extra);
       }
       if (upperBound(iColumn)) {
-        double hHat = wValue - wVec_[iColumn]*rhsU_[iColumn];
+        CoinWorkDouble hHat = wValue - wVec_[iColumn]*rhsU_[iColumn];
         value += hHat/(upperSlack_[iColumn]+extra);
       }
@@ -2667 +2741 @@
     // KKT
     for (iColumn=0;iColumn<numberTotal;iColumn++) {
-      double value = rhsC_[iColumn];
-      double zValue = rhsZ_[iColumn];
-      double wValue = rhsW_[iColumn];
+      CoinWorkDouble value = rhsC_[iColumn];
+      CoinWorkDouble zValue = rhsZ_[iColumn];
+      CoinWorkDouble wValue = rhsW_[iColumn];
       if (lowerBound(iColumn)) {
-        double gHat = zValue + zVec_[iColumn]*rhsL_[iColumn];
+        CoinWorkDouble gHat = zValue + zVec_[iColumn]*rhsL_[iColumn];
         value -= gHat/(lowerSlack_[iColumn]+extra);
       }
       if (upperBound(iColumn)) {
-        double hHat = wValue - wVec_[iColumn]*rhsU_[iColumn];
+        CoinWorkDouble hHat = wValue - wVec_[iColumn]*rhsU_[iColumn];
         value += hHat/(upperSlack_[iColumn]+extra);
       }
@@ -2684 +2758 @@
 //method: sees if looks plausible change in complementarity
 bool ClpPredictorCorrector::checkGoodMove(const bool doCorrector,
-                                          double & bestNextGap,
+                                          CoinWorkDouble & bestNextGap,
                                           bool allowIncreasingGap)
 {
-  const double beta3 = 0.99997;
+  const CoinWorkDouble beta3 = 0.99997;
   bool goodMove=false;
   int nextNumber;
   int nextNumberItems;
   int numberTotal = numberRows_+numberColumns_;
-  double returnGap=bestNextGap;
-  double nextGap=complementarityGap(nextNumber,nextNumberItems,2);
+  CoinWorkDouble returnGap=bestNextGap;
+  CoinWorkDouble nextGap=complementarityGap(nextNumber,nextNumberItems,2);
 #ifndef NO_RTTI
   ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(objective_));
@@ -2711 +2785 @@
     returnGap=nextGap;
   }
-  double step;
+  CoinWorkDouble step;
   if (actualDualStep_>actualPrimalStep_) {
     step=actualDualStep_;
@@ -2717 +2791 @@
     step=actualPrimalStep_;
   } 
-  double testValue=1.0-step*(1.0-beta3);
+  CoinWorkDouble testValue=1.0-step*(1.0-beta3);
   //testValue=0.0;
   testValue*=complementarityGap_;
@@ -2729 +2803 @@
       //step=actualPrimalStep_;
     //} 
-    double gap = bestNextGap;
+    CoinWorkDouble gap = bestNextGap;
     goodMove=checkGoodMove2(step,gap,allowIncreasingGap);
     if (goodMove)
@@ -2760 +2834 @@
     while (!goodMove) {
       pass++;
-      double gap = bestNextGap;
+      CoinWorkDouble gap = bestNextGap;
       goodMove=checkGoodMove2(step,gap,allowIncreasingGap);
       if (goodMove||pass>3) {
@@ -2793 +2867 @@
   if (goodMove) {
     //compute delta in objectives
-    double deltaObjectivePrimal=0.0;
-    double deltaObjectiveDual=
+    CoinWorkDouble deltaObjectivePrimal=0.0;
+    CoinWorkDouble deltaObjectiveDual=
       innerProduct(deltaY_,numberRows_,
                    rhsFixRegion_);
-    double error=0.0;
-    double * workArray = workArray_;
+    CoinWorkDouble error=0.0;
+    CoinWorkDouble * workArray = workArray_;
     CoinZeroN(workArray,numberColumns_);
     CoinMemcpyN(deltaY_,numberRows_,workArray+numberColumns_);
     matrix_->transposeTimes(-1.0,deltaY_,workArray);
-    //double sumPerturbCost=0.0;
+    //CoinWorkDouble sumPerturbCost=0.0;
     for (int iColumn=0;iColumn<numberTotal;iColumn++) {
       if (!flagged(iColumn)) {
@@ -2813 +2887 @@
           deltaObjectiveDual-=deltaW_[iColumn]*upper_[iColumn];
         } 
-        double change = fabs(workArray_[iColumn]-deltaZ_[iColumn]+deltaW_[iColumn]);
+        CoinWorkDouble change = CoinAbs(workArray_[iColumn]-deltaZ_[iColumn]+deltaW_[iColumn]);
         error = CoinMax(change,error);
       } 
@@ -2819 +2893 @@
     } 
     //deltaObjectivePrimal+=sumPerturbCost*linearPerturbation_;
-    double testValue;
+    CoinWorkDouble testValue;
     if (error>0.0) {
       testValue=1.0e1*CoinMax(maximumDualError_,1.0e-12)/error;
@@ -2828 +2902 @@
     if (testValue<actualDualStep_&&!quadraticObj) {
       handler_->message(CLP_BARRIER_REDUCING,messages_)
-      <<"dual"<<actualDualStep_
-      << testValue
+        <<"dual"<<static_cast<double>(actualDualStep_)
+        << static_cast<double>(testValue)
       <<CoinMessageEol;
       actualDualStep_=testValue;
@@ -2838 +2912 @@
     //check change in AX not too much
     //??? could be dropped row going infeasible
-    double ratio = 1.0e1*CoinMax(maximumRHSError_,1.0e-12)/maximumRHSChange_;
+    CoinWorkDouble ratio = 1.0e1*CoinMax(maximumRHSError_,1.0e-12)/maximumRHSChange_;
     if (ratio<actualPrimalStep_) {
       handler_->message(CLP_BARRIER_REDUCING,messages_)
-      <<"primal"<<actualPrimalStep_
-      <<ratio
+        <<"primal"<<static_cast<double>(actualPrimalStep_)
+        <<static_cast<double>(ratio)
       <<CoinMessageEol;
       if (ratio>1.0e-6) {
@@ -2857 +2931 @@
 }
 //:  checks for one step size
-bool ClpPredictorCorrector::checkGoodMove2(double move,
-                                           double & bestNextGap,
+bool ClpPredictorCorrector::checkGoodMove2(CoinWorkDouble move,
+                                           CoinWorkDouble & bestNextGap,
                                            bool allowIncreasingGap)
 {
-  double complementarityMultiplier =1.0/numberComplementarityPairs_;
-  const double gamma = 1.0e-8;
-  const double gammap = 1.0e-8;
-  double gammad = 1.0e-8;
+  CoinWorkDouble complementarityMultiplier =1.0/numberComplementarityPairs_;
+#if CLP_LONG_CHOLESKY<2
+  const CoinWorkDouble gamma = 1.0e-8;
+  const CoinWorkDouble gammap = 1.0e-8;
+  CoinWorkDouble gammad = 1.0e-8;
+#else
+  const CoinWorkDouble gamma = 1.0e-8;
+  const CoinWorkDouble gammap = 1.0e-8;
+  CoinWorkDouble gammad = 1.0e-8;
+#endif
   int nextNumber;
   int nextNumberItems;
-  double nextGap=complementarityGap(nextNumber,nextNumberItems,2);
+  CoinWorkDouble nextGap=complementarityGap(nextNumber,nextNumberItems,2);
   if (nextGap>bestNextGap&&!allowIncreasingGap)
     return false;
-  double lowerBoundGap = gamma*nextGap*complementarityMultiplier;
+  CoinWorkDouble lowerBoundGap = gamma*nextGap*complementarityMultiplier;
   bool goodMove=true;
   int iColumn;
@@ -2876 +2956 @@
     if (!flagged(iColumn)) {
       if (lowerBound(iColumn)) {
-        double part1=lowerSlack_[iColumn]+actualPrimalStep_*deltaSL_[iColumn];
-        double part2=zVec_[iColumn]+actualDualStep_*deltaZ_[iColumn];
+        CoinWorkDouble part1=lowerSlack_[iColumn]+actualPrimalStep_*deltaSL_[iColumn];
+        CoinWorkDouble part2=zVec_[iColumn]+actualDualStep_*deltaZ_[iColumn];
         if (part1*part2<lowerBoundGap) {
           goodMove=false;
@@ -2884 +2964 @@
       } 
       if (upperBound(iColumn)) {
-        double part1=upperSlack_[iColumn]+actualPrimalStep_*deltaSU_[iColumn];
-        double part2=wVec_[iColumn]+actualDualStep_*deltaW_[iColumn];
+        CoinWorkDouble part1=upperSlack_[iColumn]+actualPrimalStep_*deltaSU_[iColumn];
+        CoinWorkDouble part2=wVec_[iColumn]+actualDualStep_*deltaW_[iColumn];
         if (part1*part2<lowerBoundGap) {
           goodMove=false;
@@ -2893 +2973 @@
     } 
   } 
-   double * nextDj=NULL;
-   double maximumDualError = maximumDualError_;
+   CoinWorkDouble * nextDj=NULL;
+   CoinWorkDouble maximumDualError = maximumDualError_;
 #ifndef NO_RTTI
   ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(objective_));
@@ -2902 +2982 @@
     quadraticObj = (static_cast< ClpQuadraticObjective*>(objective_));
 #endif
+  CoinWorkDouble * dualArray = reinterpret_cast<CoinWorkDouble *>(dual_);
   if (quadraticObj) {
     // change gammad
     gammad=1.0e-4;
-    double gamma2 = gamma_*gamma_;
-    nextDj = new double [numberColumns_];
-    double * nextSolution = new double [numberColumns_];
+    CoinWorkDouble gamma2 = gamma_*gamma_;
+    nextDj = new CoinWorkDouble [numberColumns_];
+    CoinWorkDouble * nextSolution = new CoinWorkDouble [numberColumns_];
     // put next primal into nextSolution
     for ( iColumn=0;iColumn<numberColumns_;iColumn++) {
@@ -2919 +3000 @@
     // do reduced costs
     CoinMemcpyN(cost_,numberColumns_,nextDj);
-    matrix_->transposeTimes(-1.0,dual_,nextDj);
+    matrix_->transposeTimes(-1.0,dualArray,nextDj);
     matrix_->transposeTimes(-actualDualStep_,deltaY_,nextDj);
     quadraticDjs(nextDj,nextSolution,1.0);
@@ -2927 +3008 @@
     for (int iColumn=0;iColumn<numberColumns_;iColumn++) {
       if (!fixedOrFree(iColumn)) {
-        double newZ=0.0;
-        double newW=0.0;
+        CoinWorkDouble newZ=0.0;
+        CoinWorkDouble newW=0.0;
         if (lowerBound(iColumn)) {
           newZ=zVec_[iColumn]+actualDualStep_*deltaZ_[iColumn];
@@ -2936 +3017 @@
         } 
         if (columnQuadraticLength[iColumn]) {
-          double gammaTerm = gamma2;
+          CoinWorkDouble gammaTerm = gamma2;
           if (primalR_)
             gammaTerm += primalR_[iColumn];
-          //double dualInfeasibility=
+          //CoinWorkDouble dualInfeasibility=
           //dj_[iColumn]-zVec_[iColumn]+wVec_[iColumn]
           //+gammaTerm*solution_[iColumn];
-          double newInfeasibility=
+          CoinWorkDouble newInfeasibility=
             nextDj[iColumn]-newZ+newW
             +gammaTerm*(solution_[iColumn]+actualPrimalStep_*deltaX_[iColumn]);
           maximumDualError = CoinMax(maximumDualError,newInfeasibility);
-          //if (fabs(newInfeasibility)>CoinMax(2000.0*maximumDualError_,1.0e-2)) {
+          //if (CoinAbs(newInfeasibility)>CoinMax(2000.0*maximumDualError_,1.0e-2)) {
           //if (dualInfeasibility*newInfeasibility<0.0) {
           //  printf("%d current %g next %g\n",iColumn,dualInfeasibility,
@@ -2962 +3043 @@
     solutionNorm_=rhsNorm_;
   } 
-  double errorCheck=maximumRHSError_/solutionNorm_;
+  CoinWorkDouble errorCheck=maximumRHSError_/solutionNorm_;
   if (errorCheck<maximumBoundInfeasibility_) {
     errorCheck=maximumBoundInfeasibility_;
@@ -2987 +3068 @@
 // updateSolution.  Updates solution at end of iteration
 //returns number fixed
-int ClpPredictorCorrector::updateSolution(double nextGap)
+int ClpPredictorCorrector::updateSolution(CoinWorkDouble nextGap)
 {
+  CoinWorkDouble * dualArray = reinterpret_cast<CoinWorkDouble *>(dual_);
   int numberTotal = numberRows_+numberColumns_;
   //update pi
-  multiplyAdd(deltaY_,numberRows_,actualDualStep_,dual_,1.0);
+  multiplyAdd(deltaY_,numberRows_,actualDualStep_,dualArray,1.0);
   CoinZeroN(errorRegion_,numberRows_);
   CoinZeroN(rhsFixRegion_,numberRows_);
-  double maximumRhsInfeasibility=0.0;
-  double maximumBoundInfeasibility=0.0;
-  double maximumDualError=1.0e-12;
-  double primalObjectiveValue=0.0;
-  double dualObjectiveValue=0.0;
-  double solutionNorm=1.0e-12;
+  CoinWorkDouble maximumRhsInfeasibility=0.0;
+  CoinWorkDouble maximumBoundInfeasibility=0.0;
+  CoinWorkDouble maximumDualError=1.0e-12;
+  CoinWorkDouble primalObjectiveValue=0.0;
+  CoinWorkDouble dualObjectiveValue=0.0;
+  CoinWorkDouble solutionNorm=1.0e-12;
   int numberKilled=0;
-  double freeMultiplier=1.0e6;
-  double trueNorm =diagonalNorm_/diagonalScaleFactor_;
+  CoinWorkDouble freeMultiplier=1.0e6;
+  CoinWorkDouble trueNorm =diagonalNorm_/diagonalScaleFactor_;
   if (freeMultiplier<trueNorm) {
     freeMultiplier=trueNorm;
@@ -3010 +3092 @@
   } 
   freeMultiplier=0.5/freeMultiplier;
-  double condition = fabs(cholesky_->choleskyCondition());
+  CoinWorkDouble condition = CoinAbs(cholesky_->choleskyCondition());
   bool caution;
   if ((condition<1.0e10&&trueNorm<1.0e12)||numberIterations_<20) {
@@ -3017 +3099 @@
     caution=true;
   } 
-  double extra=eExtra;
-  const double largeFactor=1.0e2;
-  double largeGap=largeFactor*solutionNorm_;
+  CoinWorkDouble extra=eExtra;
+  const CoinWorkDouble largeFactor=1.0e2;
+  CoinWorkDouble largeGap=largeFactor*solutionNorm_;
   if (largeGap<largeFactor) {
     largeGap=largeFactor;
   } 
-  double dualFake=0.0;
-  double dualTolerance =  dblParam_[ClpDualTolerance];
+  CoinWorkDouble dualFake=0.0;
+  CoinWorkDouble dualTolerance =  dblParam_[ClpDualTolerance];
   dualTolerance=dualTolerance/scaleFactor_;
   if (dualTolerance<1.0e-12) {
     dualTolerance=1.0e-12;
   } 
-  double offsetObjective=0.0;
-  const double killTolerance=primalTolerance();
-  double qDiagonal;
+  CoinWorkDouble offsetObjective=0.0;
+  const CoinWorkDouble killTolerance=primalTolerance();
+  CoinWorkDouble qDiagonal;
+#if CLP_LONG_CHOLESKY<2
   if (mu_<1.0) {
     qDiagona