- Timestamp:
- May 23, 2009 3:47:00 AM (12 years ago)
- Location:
- trunk/Clp/src
- Files:
-
- 22 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/Clp/src/CbcOrClpParam.cpp
r1365 r1367 2868 2868 ); 2869 2869 parameters[numberParameters-1].setIntValue(1); 2870 #ifdef CBC_KEEP_DEPRECATED 2870 2871 parameters[numberParameters++]= 2871 2872 CbcOrClpParam("strengthen","Create strengthened problem", … … 2875 2876 "This creates a new problem by applying the root node cuts. All tight constraints \ 2876 2877 will be in resulting problem" 2877 ); 2878 ); 2879 #endif 2878 2880 parameters[numberParameters++]= 2879 2881 CbcOrClpParam("strong!Branching","Number of variables to look at in strong branching", -
trunk/Clp/src/ClpCholeskyBase.cpp
r1321 r1367 86 86 #else 87 87 // actually long double 88 workDouble_ = (double *) ClpCopyOfArray((longWork *) rhs.workDouble_,numberRows_);88 workDouble_ = reinterpret_cast<double *> (ClpCopyOfArray(reinterpret_cast<CoinWorkDouble *> (rhs.workDouble_),numberRows_)); 89 89 #endif 90 90 link_ = ClpCopyOfArray(rhs.link_,numberRows_); … … 171 171 #else 172 172 // actually long double 173 workDouble_ = (double *) ClpCopyOfArray((longWork *) rhs.workDouble_,numberRows_);173 workDouble_ = reinterpret_cast<double *> (ClpCopyOfArray(reinterpret_cast<CoinWorkDouble *> (rhs.workDouble_),numberRows_)); 174 174 #endif 175 175 link_ = ClpCopyOfArray(rhs.link_,numberRows_); … … 195 195 region1 is rows+columns, region2 is rows */ 196 196 void 197 ClpCholeskyBase::solveKKT ( double * region1, double * region2, const double * diagonal,198 double diagonalScaleFactor)197 ClpCholeskyBase::solveKKT (CoinWorkDouble * region1, CoinWorkDouble * region2, const CoinWorkDouble * diagonal, 198 CoinWorkDouble diagonalScaleFactor) 199 199 { 200 200 if (!doKKT_) { … … 202 202 int numberColumns = model_->numberColumns(); 203 203 int numberTotal = numberRows_+numberColumns; 204 double * region1Save = new double[numberTotal];204 CoinWorkDouble * region1Save = new CoinWorkDouble[numberTotal]; 205 205 for (iColumn=0;iColumn<numberTotal;iColumn++) { 206 206 region1[iColumn] *= diagonal[iColumn]; … … 209 209 multiplyAdd(region1+numberColumns,numberRows_,-1.0,region2,1.0); 210 210 model_->clpMatrix()->times(1.0,region1,region2); 211 double maximumRHS = maximumAbsElement(region2,numberRows_);212 double scale=1.0;213 double unscale=1.0;211 CoinWorkDouble maximumRHS = maximumAbsElement(region2,numberRows_); 212 CoinWorkDouble scale=1.0; 213 CoinWorkDouble unscale=1.0; 214 214 if (maximumRHS>1.0e-30) { 215 215 if (maximumRHS<=0.5) { 216 double factor=2.0;216 CoinWorkDouble factor=2.0; 217 217 while (maximumRHS<=0.5) { 218 218 maximumRHS*=factor; … … 220 220 } /* endwhile */ 221 221 } else if (maximumRHS>=2.0&&maximumRHS<=COIN_DBL_MAX) { 222 double factor=0.5;222 CoinWorkDouble factor=0.5; 223 223 while (maximumRHS>=2.0) { 224 224 maximumRHS*=factor; … … 246 246 int numberColumns = model_->numberColumns(); 247 247 int numberTotal = numberColumns + numberRowsModel; 248 double * array = new double [numberRows_];248 CoinWorkDouble * array = new CoinWorkDouble [numberRows_]; 249 249 CoinMemcpyN(region1,numberTotal,array); 250 250 CoinMemcpyN(region2,numberRowsModel,array+numberTotal); … … 253 253 int iRow; 254 254 for (iRow=0;iRow<numberTotal;iRow++) { 255 if (rowsDropped_[iRow]&& fabs(array[iRow])>1.0e-8) {255 if (rowsDropped_[iRow]&&CoinAbs(array[iRow])>1.0e-8) { 256 256 printf("row region1 %d dropped %g\n",iRow,array[iRow]); 257 257 } 258 258 } 259 259 for (;iRow<numberRows_;iRow++) { 260 if (rowsDropped_[iRow]&& fabs(array[iRow])>1.0e-8) {260 if (rowsDropped_[iRow]&&CoinAbs(array[iRow])>1.0e-8) { 261 261 printf("row region2 %d dropped %g\n",iRow,array[iRow]); 262 262 } … … 1237 1237 #else 1238 1238 // actually long double 1239 workDouble_ = (double *) (new longWork[numberRows_]);1239 workDouble_ = reinterpret_cast<double *> (new CoinWorkDouble[numberRows_]); 1240 1240 #endif 1241 1241 diagonal_ = new longDouble[numberRows_]; … … 1521 1521 /* Factorize - filling in rowsDropped and returning number dropped */ 1522 1522 int 1523 ClpCholeskyBase::factorize(const double * diagonal, int * rowsDropped)1523 ClpCholeskyBase::factorize(const CoinWorkDouble * diagonal, int * rowsDropped) 1524 1524 { 1525 1525 const CoinBigIndex * columnStart = model_->clpMatrix()->getVectorStarts(); … … 1533 1533 int numberColumns=model_->clpMatrix()->getNumCols(); 1534 1534 //perturbation 1535 longDouble perturbation=model_->diagonalPerturbation()*model_->diagonalNorm();1535 CoinWorkDouble perturbation=model_->diagonalPerturbation()*model_->diagonalNorm(); 1536 1536 //perturbation=perturbation*perturbation*100000000.0; 1537 1537 if (perturbation>1.0) { … … 1540 1540 std::cout <<"large perturbation "<<perturbation<<std::endl; 1541 1541 #endif 1542 perturbation= sqrt(perturbation);;1542 perturbation=CoinSqrt(perturbation); 1543 1543 perturbation=1.0; 1544 1544 } … … 1548 1548 CoinZeroN(work,numberRows_); 1549 1549 int newDropped=0; 1550 double largest=1.0;1551 double smallest=COIN_DBL_MAX;1550 CoinWorkDouble largest=1.0; 1551 CoinWorkDouble smallest=COIN_DBL_MAX; 1552 1552 int numberDense=0; 1553 1553 if (!doKKT_) { 1554 const double * diagonalSlack = diagonal + numberColumns;1554 const CoinWorkDouble * diagonalSlack = diagonal + numberColumns; 1555 1555 if (dense_) 1556 1556 numberDense=dense_->numberRows(); … … 1577 1577 } 1578 1578 } 1579 longDouble delta2 = model_->delta(); // add delta*delta to diagonal1579 CoinWorkDouble delta2 = model_->delta(); // add delta*delta to diagonal 1580 1580 delta2 *= delta2; 1581 1581 // largest in initial matrix 1582 double largest2=1.0e-20;1582 CoinWorkDouble largest2=1.0e-20; 1583 1583 for (iRow=0;iRow<numberRows_;iRow++) { 1584 1584 longDouble * put = sparseFactor_+choleskyStart_[iRow]; … … 1597 1597 CoinBigIndex start=columnStart[iColumn]; 1598 1598 CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn]; 1599 longDouble multiplier = diagonal[iColumn]*elementByRow[k];1599 CoinWorkDouble multiplier = diagonal[iColumn]*elementByRow[k]; 1600 1600 for (CoinBigIndex j=start;j<end;j++) { 1601 1601 int jRow=row[j]; 1602 1602 int jNewRow = permuteInverse_[jRow]; 1603 1603 if (jNewRow>=iRow&&!rowsDropped_[jRow]) { 1604 longDouble value=element[j]*multiplier;1604 CoinWorkDouble value=element[j]*multiplier; 1605 1605 work[jNewRow] += value; 1606 1606 } … … 1609 1609 } 1610 1610 diagonal_[iRow]=work[iRow]; 1611 largest2 = CoinMax(largest2, fabs(work[iRow]));1611 largest2 = CoinMax(largest2,CoinAbs(work[iRow])); 1612 1612 work[iRow]=0.0; 1613 1613 int j; … … 1615 1615 int jRow = which[j]; 1616 1616 put[j]=work[jRow]; 1617 largest2 = CoinMax(largest2, fabs(work[jRow]));1617 largest2 = CoinMax(largest2,CoinAbs(work[jRow])); 1618 1618 work[jRow]=0.0; 1619 1619 } … … 1629 1629 //check sizes 1630 1630 largest2*=1.0e-20; 1631 largest = CoinMin(largest2, 1.0e-11);1631 largest = CoinMin(largest2,CHOL_SMALL_VALUE); 1632 1632 int numberDroppedBefore=0; 1633 1633 for (iRow=0;iRow<numberRows_;iRow++) { … … 1636 1636 rowsDropped[iRow]=dropped; 1637 1637 if (!dropped) { 1638 longDouble diagonal = diagonal_[iRow];1638 CoinWorkDouble diagonal = diagonal_[iRow]; 1639 1639 if (diagonal>largest2) { 1640 1640 diagonal_[iRow]=diagonal+perturbation; … … 1663 1663 for ( i=0;i<numberRows_;i++) { 1664 1664 assert (diagonal_[i]>=0.0); 1665 diagonal_[i]= sqrt(diagonal_[i]);1665 diagonal_[i]=CoinSqrt(diagonal_[i]); 1666 1666 } 1667 1667 // Update dense columns (just L) … … 1673 1673 a[j]=0.0; 1674 1674 } 1675 solveLong(a,1); 1675 for (i=0;i<numberRows_;i++) { 1676 int iRow = permute_[i]; 1677 workDouble_[i] = a[iRow]; 1678 } 1679 for (i=0;i<numberRows_;i++) { 1680 CoinWorkDouble value=workDouble_[i]; 1681 CoinBigIndex offset = indexStart_[i]-choleskyStart_[i]; 1682 CoinBigIndex j; 1683 for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) { 1684 int iRow = choleskyRow_[j+offset]; 1685 workDouble_[iRow] -= sparseFactor_[j]*value; 1686 } 1687 } 1688 for (i=0;i<numberRows_;i++) { 1689 int iRow = permute_[i]; 1690 a[iRow]=workDouble_[i]*diagonal_[i]; 1691 } 1676 1692 } 1677 1693 dense_->resetRowsDropped(); … … 1682 1698 const longDouble * a = denseColumn_+i*numberRows_; 1683 1699 // do diagonal 1684 longDouble value = denseDiagonal[i];1700 CoinWorkDouble value = denseDiagonal[i]; 1685 1701 const longDouble * b = denseColumn_+i*numberRows_; 1686 1702 for (int k=0;k<numberRows_;k++) … … 1688 1704 denseDiagonal[i]=value; 1689 1705 for (int j=i+1;j<numberDense;j++) { 1690 longDouble value = 0.0;1706 CoinWorkDouble value = 0.0; 1691 1707 const longDouble * b = denseColumn_+j*numberRows_; 1692 1708 for (int k=0;k<numberRows_;k++) … … 1778 1794 } 1779 1795 if (permuted) { 1780 longDouble delta2 = model_->delta(); // add delta*delta to bottom1796 CoinWorkDouble delta2 = model_->delta(); // add delta*delta to bottom 1781 1797 delta2 *= delta2; 1782 1798 // Need to permute - ugly … … 1788 1804 if (iOriginalRow<numberColumns) { 1789 1805 iColumn=iOriginalRow; 1790 longDouble value = diagonal[iColumn];1791 if ( fabs(value)>1.0e-100) {1806 CoinWorkDouble value = diagonal[iColumn]; 1807 if (CoinAbs(value)>1.0e-100) { 1792 1808 value = 1.0/value; 1793 largest = CoinMax(largest, fabs(value));1809 largest = CoinMax(largest,CoinAbs(value)); 1794 1810 diagonal_[iRow] = -value; 1795 1811 CoinBigIndex start=columnStart[iColumn]; … … 1800 1816 if (kRow>iRow) { 1801 1817 work[kRow]=element[j]; 1802 largest = CoinMax(largest, fabs(element[j]));1818 largest = CoinMax(largest,CoinAbs(element[j])); 1803 1819 } 1804 1820 } … … 1807 1823 } 1808 1824 } else if (iOriginalRow<numberTotal) { 1809 longDouble value = diagonal[iOriginalRow];1810 if ( fabs(value)>1.0e-100) {1825 CoinWorkDouble value = diagonal[iOriginalRow]; 1826 if (CoinAbs(value)>1.0e-100) { 1811 1827 value = 1.0/value; 1812 largest = CoinMax(largest, fabs(value));1828 largest = CoinMax(largest,CoinAbs(value)); 1813 1829 } else { 1814 1830 value = 1.0e100; … … 1828 1844 if (jNewRow>iRow) { 1829 1845 work[jNewRow]=elementByRow[j]; 1830 largest = CoinMax(largest, fabs(elementByRow[j]));1846 largest = CoinMax(largest,CoinAbs(elementByRow[j])); 1831 1847 } 1832 1848 } … … 1857 1873 CoinBigIndex j; 1858 1874 iColumn=iOriginalRow; 1859 longDouble value = diagonal[iColumn];1860 if ( fabs(value)>1.0e-100) {1875 CoinWorkDouble value = diagonal[iColumn]; 1876 if (CoinAbs(value)>1.0e-100) { 1861 1877 value = 1.0/value; 1862 1878 for (j=columnQuadraticStart[iColumn]; … … 1870 1886 } 1871 1887 } 1872 largest = CoinMax(largest, fabs(value));1888 largest = CoinMax(largest,CoinAbs(value)); 1873 1889 diagonal_[iRow] = -value; 1874 1890 CoinBigIndex start=columnStart[iColumn]; … … 1879 1895 if (kRow>iRow) { 1880 1896 work[kRow]=element[j]; 1881 largest = CoinMax(largest, fabs(element[j]));1897 largest = CoinMax(largest,CoinAbs(element[j])); 1882 1898 } 1883 1899 } … … 1886 1902 } 1887 1903 } else if (iOriginalRow<numberTotal) { 1888 longDouble value = diagonal[iOriginalRow];1889 if ( fabs(value)>1.0e-100) {1904 CoinWorkDouble value = diagonal[iOriginalRow]; 1905 if (CoinAbs(value)>1.0e-100) { 1890 1906 value = 1.0/value; 1891 largest = CoinMax(largest, fabs(value));1907 largest = CoinMax(largest,CoinAbs(value)); 1892 1908 } else { 1893 1909 value = 1.0e100; … … 1907 1923 if (jNewRow>iRow) { 1908 1924 work[jNewRow]=elementByRow[j]; 1909 largest = CoinMax(largest, fabs(elementByRow[j]));1925 largest = CoinMax(largest,CoinAbs(elementByRow[j])); 1910 1926 } 1911 1927 } … … 1931 1947 longDouble * put = sparseFactor_+choleskyStart_[iColumn]; 1932 1948 int * which = choleskyRow_+indexStart_[iColumn]; 1933 longDouble value = diagonal[iColumn];1934 if ( fabs(value)>1.0e-100) {1949 CoinWorkDouble value = diagonal[iColumn]; 1950 if (CoinAbs(value)>1.0e-100) { 1935 1951 value = 1.0/value; 1936 largest = CoinMax(largest, fabs(value));1952 largest = CoinMax(largest,CoinAbs(value)); 1937 1953 diagonal_[iColumn] = -value; 1938 1954 CoinBigIndex start=columnStart[iColumn]; … … 1942 1958 //sparseFactor_[numberElements++]=element[j]; 1943 1959 work[row[j]+numberTotal]=element[j]; 1944 largest = CoinMax(largest, fabs(element[j]));1960 largest = CoinMax(largest,CoinAbs(element[j])); 1945 1961 } 1946 1962 } else { … … 1965 1981 int * which = choleskyRow_+indexStart_[iColumn]; 1966 1982 int number = choleskyStart_[iColumn+1]-choleskyStart_[iColumn]; 1967 longDouble value = diagonal[iColumn];1983 CoinWorkDouble value = diagonal[iColumn]; 1968 1984 CoinBigIndex j; 1969 if ( fabs(value)>1.0e-100) {1985 if (CoinAbs(value)>1.0e-100) { 1970 1986 value = 1.0/value; 1971 1987 for (j=columnQuadraticStart[iColumn]; … … 1978 1994 } 1979 1995 } 1980 largest = CoinMax(largest, fabs(value));1996 largest = CoinMax(largest,CoinAbs(value)); 1981 1997 diagonal_[iColumn] = -value; 1982 1998 CoinBigIndex start=columnStart[iColumn]; … … 1984 2000 for (j=start;j<end;j++) { 1985 2001 work[row[j]+numberTotal]=element[j]; 1986 largest = CoinMax(largest, fabs(element[j]));2002 largest = CoinMax(largest,CoinAbs(element[j])); 1987 2003 } 1988 2004 for (j=0;j<number;j++) { … … 2005 2021 longDouble * put = sparseFactor_+choleskyStart_[iColumn]; 2006 2022 int * which = choleskyRow_+indexStart_[iColumn]; 2007 longDouble value = diagonal[iColumn];2008 if ( fabs(value)>1.0e-100) {2023 CoinWorkDouble value = diagonal[iColumn]; 2024 if (CoinAbs(value)>1.0e-100) { 2009 2025 value = 1.0/value; 2010 largest = CoinMax(largest, fabs(value));2026 largest = CoinMax(largest,CoinAbs(value)); 2011 2027 } else { 2012 2028 value = 1.0e100; … … 2023 2039 } 2024 2040 // Finish diagonal 2025 longDouble delta2 = model_->delta(); // add delta*delta to bottom2041 CoinWorkDouble delta2 = model_->delta(); // add delta*delta to bottom 2026 2042 delta2 *= delta2; 2027 2043 for (iRow=0;iRow<numberRowsModel;iRow++) { … … 2037 2053 //check sizes 2038 2054 largest*=1.0e-20; 2039 largest = CoinMin(largest, 1.0e-11);2055 largest = CoinMin(largest,CHOL_SMALL_VALUE); 2040 2056 doubleParameters_[10]=CoinMax(1.0e-20,largest); 2041 2057 integerParameters_[20]=0; … … 2056 2072 // Should save adjustments in ..R_ 2057 2073 int n1=0,n2=0; 2058 double * primalR = model_->primalR();2059 double * dualR = model_->dualR();2074 CoinWorkDouble * primalR = model_->primalR(); 2075 CoinWorkDouble * dualR = model_->dualR(); 2060 2076 for (iRow=0;iRow<numberTotal;iRow++) { 2061 2077 if (rowsDropped2[iRow]) { … … 2093 2109 ClpCholeskyBase::factorizePart2(int * rowsDropped) 2094 2110 { 2095 longDouble largest=doubleParameters_[3];2096 longDouble smallest=doubleParameters_[4];2111 CoinWorkDouble largest=doubleParameters_[3]; 2112 CoinWorkDouble smallest=doubleParameters_[4]; 2097 2113 int numberDropped=integerParameters_[20]; 2098 2114 // probably done before … … 2159 2175 for ( jRow=lastRow;jRow<iRow;jRow++) { 2160 2176 int jCount = jRow-lastRow; 2161 longDouble diagonalValue = diagonal_[jRow];2177 CoinWorkDouble diagonalValue = diagonal_[jRow]; 2162 2178 CoinBigIndex start=choleskyStart_[jRow]; 2163 2179 CoinBigIndex end=choleskyStart_[jRow+1]; … … 2165 2181 jCount--; 2166 2182 CoinBigIndex get = choleskyStart_[kRow]+jCount; 2167 longDouble a_jk = sparseFactor_[get];2168 longDouble value1 = d[kRow]*a_jk;2183 CoinWorkDouble a_jk = sparseFactor_[get]; 2184 CoinWorkDouble value1 = d[kRow]*a_jk; 2169 2185 diagonalValue -= a_jk*value1; 2170 2186 for (CoinBigIndex j=start;j<end;j++) … … 2222 2238 } 2223 2239 // for each column L[*,kRow] that affects L[*,iRow] 2224 longDouble diagonalValue=diagonal_[iRow];2240 CoinWorkDouble diagonalValue=diagonal_[iRow]; 2225 2241 int nextRow = link_[iRow]; 2226 2242 int kRow=0; … … 2234 2250 CoinBigIndex end=choleskyStart_[kRow+1]; 2235 2251 assert(k<end); 2236 longDouble a_ik=sparseFactor_[k++];2237 longDouble value1=d[kRow]*a_ik;2252 CoinWorkDouble a_ik=sparseFactor_[k++]; 2253 CoinWorkDouble value1=d[kRow]*a_ik; 2238 2254 // update first 2239 2255 first[kRow]=k; … … 2259 2275 CoinBigIndex j=first[kkRow]; 2260 2276 //int iiRow = choleskyRow_[j+indexStart_[kkRow]-choleskyStart_[kkRow]]; 2261 longDouble a = sparseFactor_[j];2262 longDouble dValue = d[kkRow]*a;2277 CoinWorkDouble a = sparseFactor_[j]; 2278 CoinWorkDouble dValue = d[kkRow]*a; 2263 2279 diagonalValue -= a*dValue; 2264 2280 work[kkRow]=dValue; … … 2271 2287 for (int i=0;i<length;i++) { 2272 2288 int lRow = choleskyRow_[currentIndex++]; 2273 longDouble t0 = work[lRow];2289 CoinWorkDouble t0 = work[lRow]; 2274 2290 for (int kkRow=kRow;kkRow<last;kkRow++) { 2275 2291 CoinBigIndex j = first[kkRow]+i; … … 2340 2356 for (CoinBigIndex j=start;j<end;j++) { 2341 2357 int jRow = choleskyRow_[j+offset]; 2342 longDouble value = sparseFactor_[j]-work[jRow];2358 CoinWorkDouble value = sparseFactor_[j]-work[jRow]; 2343 2359 work[jRow]=0.0; 2344 2360 sparseFactor_[j]= diagonalValue*value; … … 2395 2411 CoinBigIndex offset = indexStart_[iRow]-choleskyStart_[iRow]; 2396 2412 if (clique_[iRow]<2) { 2397 longDouble dValue=d[iRow];2413 CoinWorkDouble dValue=d[iRow]; 2398 2414 for (CoinBigIndex k=start;k<end;k++) { 2399 2415 int kRow = choleskyRow_[k+offset]; 2400 2416 assert(kRow>=firstDense_); 2401 longDouble a_ik=sparseFactor_[k];2402 longDouble value1=dValue*a_ik;2417 CoinWorkDouble a_ik=sparseFactor_[k]; 2418 CoinWorkDouble value1=dValue*a_ik; 2403 2419 diagonal_[kRow] -= value1*a_ik; 2404 2420 CoinBigIndex base = choleskyStart_[kRow]-kRow-1; 2405 2421 for (CoinBigIndex j=k+1;j<end;j++) { 2406 2422 int jRow=choleskyRow_[j+offset]; 2407 longDouble a_jk = sparseFactor_[j];2423 CoinWorkDouble a_jk = sparseFactor_[j]; 2408 2424 sparseFactor_[base+jRow] -= a_jk*value1; 2409 2425 } … … 2411 2427 } else if (clique_[iRow]<3) { 2412 2428 // do as pair 2413 longDouble dValue0=d[iRow];2414 longDouble dValue1=d[iRow+1];2429 CoinWorkDouble dValue0=d[iRow]; 2430 CoinWorkDouble dValue1=d[iRow+1]; 2415 2431 int offset1 = first[iRow+1]-start; 2416 2432 // skip row … … 2419 2435 int kRow = choleskyRow_[k+offset]; 2420 2436 assert(kRow>=firstDense_); 2421 longDouble a_ik0=sparseFactor_[k];2422 longDouble value0=dValue0*a_ik0;2423 longDouble a_ik1=sparseFactor_[k+offset1];2424 longDouble value1=dValue1*a_ik1;2437 CoinWorkDouble a_ik0=sparseFactor_[k]; 2438 CoinWorkDouble value0=dValue0*a_ik0; 2439 CoinWorkDouble a_ik1=sparseFactor_[k+offset1]; 2440 CoinWorkDouble value1=dValue1*a_ik1; 2425 2441 diagonal_[kRow] -= value0*a_ik0 + value1*a_ik1; 2426 2442 CoinBigIndex base = choleskyStart_[kRow]-kRow-1; 2427 2443 for (CoinBigIndex j=k+1;j<end;j++) { 2428 2444 int jRow=choleskyRow_[j+offset]; 2429 longDouble a_jk0 = sparseFactor_[j];2430 longDouble a_jk1 = sparseFactor_[j+offset1];2445 CoinWorkDouble a_jk0 = sparseFactor_[j]; 2446 CoinWorkDouble a_jk1 = sparseFactor_[j+offset1]; 2431 2447 sparseFactor_[base+jRow] -= a_jk0*value0+a_jk1*value1; 2432 2448 } … … 2441 2457 // maybe later get fancy on big cliques and do transpose copy 2442 2458 // seems only worth going to 3 on Intel 2443 longDouble dValue0=d[iRow];2444 longDouble dValue1=d[iRow+1];2445 longDouble dValue2=d[iRow+2];2459 CoinWorkDouble dValue0=d[iRow]; 2460 CoinWorkDouble dValue1=d[iRow+1]; 2461 CoinWorkDouble dValue2=d[iRow+2]; 2446 2462 // get offsets and skip rows 2447 2463 int offset1 = first[++iRow]-start; … … 2450 2466 int kRow = choleskyRow_[k+offset]; 2451 2467 assert(kRow>=firstDense_); 2452 double diagonalValue=diagonal_[kRow];2453 longDouble a_ik0=sparseFactor_[k];2454 longDouble value0=dValue0*a_ik0;2455 longDouble a_ik1=sparseFactor_[k+offset1];2456 longDouble value1=dValue1*a_ik1;2457 longDouble a_ik2=sparseFactor_[k+offset2];2458 longDouble value2=dValue2*a_ik2;2468 CoinWorkDouble diagonalValue=diagonal_[kRow]; 2469 CoinWorkDouble a_ik0=sparseFactor_[k]; 2470 CoinWorkDouble value0=dValue0*a_ik0; 2471 CoinWorkDouble a_ik1=sparseFactor_[k+offset1]; 2472 CoinWorkDouble value1=dValue1*a_ik1; 2473 CoinWorkDouble a_ik2=sparseFactor_[k+offset2]; 2474 CoinWorkDouble value2=dValue2*a_ik2; 2459 2475 CoinBigIndex base = choleskyStart_[kRow]-kRow-1; 2460 2476 diagonal_[kRow] = diagonalValue - value0*a_ik0 - value1*a_ik1 - value2*a_ik2; 2461 2477 for (CoinBigIndex j=k+1;j<end;j++) { 2462 2478 int jRow=choleskyRow_[j+offset]; 2463 longDouble a_jk0 = sparseFactor_[j];2464 longDouble a_jk1 = sparseFactor_[j+offset1];2465 longDouble a_jk2 = sparseFactor_[j+offset2];2479 CoinWorkDouble a_jk0 = sparseFactor_[j]; 2480 CoinWorkDouble a_jk1 = sparseFactor_[j+offset1]; 2481 CoinWorkDouble a_jk2 = sparseFactor_[j+offset2]; 2466 2482 sparseFactor_[base+jRow] -= a_jk0*value0+a_jk1*value1+a_jk2*value2; 2467 2483 } … … 2472 2488 // maybe later get fancy on big cliques and do transpose copy 2473 2489 // maybe only worth going to 3 on Intel (but may have hidden registers) 2474 longDouble dValue0=d[iRow];2475 longDouble dValue1=d[iRow+1];2476 longDouble dValue2=d[iRow+2];2477 longDouble dValue3=d[iRow+3];2490 CoinWorkDouble dValue0=d[iRow]; 2491 CoinWorkDouble dValue1=d[iRow+1]; 2492 CoinWorkDouble dValue2=d[iRow+2]; 2493 CoinWorkDouble dValue3=d[iRow+3]; 2478 2494 // get offsets and skip rows 2479 2495 int offset1 = first[++iRow]-start; … … 2483 2499 int kRow = choleskyRow_[k+offset]; 2484 2500 assert(kRow>=firstDense_); 2485 double diagonalValue=diagonal_[kRow];2486 longDouble a_ik0=sparseFactor_[k];2487 longDouble value0=dValue0*a_ik0;2488 longDouble a_ik1=sparseFactor_[k+offset1];2489 longDouble value1=dValue1*a_ik1;2490 longDouble a_ik2=sparseFactor_[k+offset2];2491 longDouble value2=dValue2*a_ik2;2492 longDouble a_ik3=sparseFactor_[k+offset3];2493 longDouble value3=dValue3*a_ik3;2501 CoinWorkDouble diagonalValue=diagonal_[kRow]; 2502 CoinWorkDouble a_ik0=sparseFactor_[k]; 2503 CoinWorkDouble value0=dValue0*a_ik0; 2504 CoinWorkDouble a_ik1=sparseFactor_[k+offset1]; 2505 CoinWorkDouble value1=dValue1*a_ik1; 2506 CoinWorkDouble a_ik2=sparseFactor_[k+offset2]; 2507 CoinWorkDouble value2=dValue2*a_ik2; 2508 CoinWorkDouble a_ik3=sparseFactor_[k+offset3]; 2509 CoinWorkDouble value3=dValue3*a_ik3; 2494 2510 CoinBigIndex base = choleskyStart_[kRow]-kRow-1; 2495 2511 diagonal_[kRow] = diagonalValue - (value0*a_ik0 + value1*a_ik1 + value2*a_ik2+value3*a_ik3); 2496 2512 for (CoinBigIndex j=k+1;j<end;j++) { 2497 2513 int jRow=choleskyRow_[j+offset]; 2498 longDouble a_jk0 = sparseFactor_[j];2499 longDouble a_jk1 = sparseFactor_[j+offset1];2500 longDouble a_jk2 = sparseFactor_[j+offset2];2501 longDouble a_jk3 = sparseFactor_[j+offset3];2514 CoinWorkDouble a_jk0 = sparseFactor_[j]; 2515 CoinWorkDouble a_jk1 = sparseFactor_[j+offset1]; 2516 CoinWorkDouble a_jk2 = sparseFactor_[j+offset2]; 2517 CoinWorkDouble a_jk3 = sparseFactor_[j+offset3]; 2502 2518 sparseFactor_[base+jRow] -= a_jk0*value0+a_jk1*value1+a_jk2*value2+a_jk3*value3; 2503 2519 } … … 2520 2536 // do change; 2521 2537 int numberDense = dense_->numberRows(); 2522 double * change = new double[numberDense];2538 CoinWorkDouble * change = new CoinWorkDouble[numberDense]; 2523 2539 for (i=0;i<numberDense;i++) { 2524 2540 const longDouble * a = denseColumn_+i*numberRows_; 2525 longDouble value =0.0;2541 CoinWorkDouble value =0.0; 2526 2542 for (int iRow=0;iRow<numberRows_;iRow++) 2527 2543 value += a[iRow]*region[iRow]; … … 2529 2545 } 2530 2546 // solve 2547 #if CLP_LONG_CHOLESKY>0 2548 dense_->solveLong(change); 2549 #else 2531 2550 dense_->solve(change); 2551 #endif 2532 2552 for (i=0;i<numberDense;i++) { 2533 2553 const longDouble * a = denseColumn_+i*numberRows_; 2534 longDouble value = change[i];2554 CoinWorkDouble value = change[i]; 2535 2555 for (int iRow=0;iRow<numberRows_;iRow++) 2536 2556 region[iRow] -= value*a[iRow]; … … 2549 2569 #ifdef CLP_DEBUG 2550 2570 double * regionX=NULL; 2551 if (sizeof( longWork)!=sizeof(double)&&type==3) {2571 if (sizeof(CoinWorkDouble)!=sizeof(double)&&type==3) { 2552 2572 regionX=ClpCopyOfArray(region,numberRows_); 2553 2573 } 2554 2574 #endif 2555 longWork * work = reinterpret_cast<longWork*> (workDouble_);2575 CoinWorkDouble * work = reinterpret_cast<CoinWorkDouble *> (workDouble_); 2556 2576 int i; 2557 2577 CoinBigIndex j; … … 2563 2583 case 1: 2564 2584 for (i=0;i<numberRows_;i++) { 2565 longDouble value=work[i];2585 CoinWorkDouble value=work[i]; 2566 2586 CoinBigIndex offset = indexStart_[i]-choleskyStart_[i]; 2567 2587 for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) { … … 2578 2598 for (i=numberRows_-1;i>=0;i--) { 2579 2599 CoinBigIndex offset = indexStart_[i]-choleskyStart_[i]; 2580 longDouble value=work[i]*diagonal_[i];2600 CoinWorkDouble value=work[i]*diagonal_[i]; 2581 2601 for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) { 2582 2602 int iRow = choleskyRow_[j+offset]; … … 2591 2611 for (i=0;i<firstDense_;i++) { 2592 2612 CoinBigIndex offset = indexStart_[i]-choleskyStart_[i]; 2593 longDouble value=work[i];2613 CoinWorkDouble value=work[i]; 2594 2614 for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) { 2595 2615 int iRow = choleskyRow_[j+offset]; … … 2609 2629 #endif 2610 2630 for (i=numberRows_-1;i>=firstDense_;i--) { 2611 longDouble value=work[i];2631 CoinWorkDouble value=work[i]; 2612 2632 int iRow = permute_[i]; 2613 2633 region[iRow]=value; … … 2616 2636 for (i=firstDense_-1;i>=0;i--) { 2617 2637 CoinBigIndex offset = indexStart_[i]-choleskyStart_[i]; 2618 longDouble value=work[i]*diagonal_[i];2638 CoinWorkDouble value=work[i]*diagonal_[i]; 2619 2639 for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) { 2620 2640 int iRow = choleskyRow_[j+offset]; … … 2634 2654 double largestO=0.0; 2635 2655 for (i=0;i<numberRows_;i++) { 2636 largestO = CoinMax(largestO, fabs(regionX[i]));2656 largestO = CoinMax(largestO,CoinAbs(regionX[i])); 2637 2657 } 2638 2658 for (i=0;i<numberRows_;i++) { … … 2642 2662 for (i=0;i<firstDense_;i++) { 2643 2663 CoinBigIndex offset = indexStart_[i]-choleskyStart_[i]; 2644 longDouble value=work[i];2664 CoinWorkDouble value=work[i]; 2645 2665 for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) { 2646 2666 int iRow = choleskyRow_[j+offset]; … … 2656 2676 dense.solveLong(work+firstDense_); 2657 2677 for (i=numberRows_-1;i>=firstDense_;i--) { 2658 longDouble value=work[i];2678 CoinWorkDouble value=work[i]; 2659 2679 int iRow = permute_[i]; 2660 2680 regionX[iRow]=value; … … 2663 2683 for (i=firstDense_-1;i>=0;i--) { 2664 2684 CoinBigIndex offset = indexStart_[i]-choleskyStart_[i]; 2665 longDouble value=work[i]*diagonal_[i];2685 CoinWorkDouble value=work[i]*diagonal_[i]; 2666 2686 for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) { 2667 2687 int iRow = choleskyRow_[j+offset]; … … 2675 2695 double largestV=0.0; 2676 2696 for (i=0;i<numberRows_;i++) { 2677 largest = CoinMax(largest, fabs(region[i]-regionX[i]));2678 largestV = CoinMax(largestV, fabs(region[i]));2697 largest = CoinMax(largest,CoinAbs(region[i]-regionX[i])); 2698 largestV = CoinMax(largestV,CoinAbs(region[i])); 2679 2699 } 2680 2700 printf("largest difference %g, largest %g, largest original %g\n", … … 2684 2704 #endif 2685 2705 } 2686 /* solve - 1 just first half, 2 just second half - 3 both. 2687 If 1 and 2 then diagonal has sqrt of inverse otherwise inverse 2688 */ 2706 #if CLP_LONG_CHOLESKY 2707 /* Uses factorization to solve. */ 2689 2708 void 2690 ClpCholeskyBase::solve Long(longDouble * region, int type)2709 ClpCholeskyBase::solve (CoinWorkDouble * region) 2691 2710 { 2711 assert (!whichDense_) ; 2712 CoinWorkDouble * work = reinterpret_cast<CoinWorkDouble *> (workDouble_); 2692 2713 int i; 2693 2714 CoinBigIndex j; 2694 2715 for (i=0;i<numberRows_;i++) { 2695 2716 int iRow = permute_[i]; 2696 workDouble_[i] = region[iRow]; 2697 } 2698 switch (type) { 2699 case 1: 2700 for (i=0;i<numberRows_;i++) { 2701 longDouble value=workDouble_[i]; 2702 CoinBigIndex offset = indexStart_[i]-choleskyStart_[i]; 2703 for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) { 2704 int iRow = choleskyRow_[j+offset]; 2705 workDouble_[iRow] -= sparseFactor_[j]*value; 2706 } 2707 } 2708 for (i=0;i<numberRows_;i++) { 2709 int iRow = permute_[i]; 2710 region[iRow]=workDouble_[i]*diagonal_[i]; 2711 } 2712 break; 2713 case 2: 2714 for (i=numberRows_-1;i>=0;i--) { 2715 CoinBigIndex offset = indexStart_[i]-choleskyStart_[i]; 2716 longDouble value=workDouble_[i]*diagonal_[i]; 2717 for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) { 2718 int iRow = choleskyRow_[j+offset]; 2719 value -= sparseFactor_[j]*workDouble_[iRow]; 2720 } 2721 workDouble_[i]=value; 2717 work[i] = region[iRow]; 2718 } 2719 for (i=0;i<firstDense_;i++) { 2720 CoinBigIndex offset = indexStart_[i]-choleskyStart_[i]; 2721 CoinWorkDouble value=work[i]; 2722 for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) { 2723 int iRow = choleskyRow_[j+offset]; 2724 work[iRow] -= sparseFactor_[j]*value; 2725 } 2726 } 2727 if (firstDense_<numberRows_) { 2728 // do dense 2729 ClpCholeskyDense dense; 2730 // just borrow space 2731 int nDense = numberRows_-firstDense_; 2732 dense.reserveSpace(this,nDense); 2733 #if CLP_LONG_CHOLESKY!=1 2734 dense.solveLong(work+firstDense_); 2735 #else 2736 dense.solveLongWork(work+firstDense_); 2737 #endif 2738 for (i=numberRows_-1;i>=firstDense_;i--) { 2739 CoinWorkDouble value=work[i]; 2722 2740 int iRow = permute_[i]; 2723 2741 region[iRow]=value; 2724 2742 } 2725 break; 2726 case 3: 2727 for (i=0;i<firstDense_;i++) { 2728 CoinBigIndex offset = indexStart_[i]-choleskyStart_[i]; 2729 longDouble value=workDouble_[i]; 2730 for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) { 2731 int iRow = choleskyRow_[j+offset]; 2732 workDouble_[iRow] -= sparseFactor_[j]*value; 2733 } 2734 } 2735 if (firstDense_<numberRows_) { 2736 // do dense 2737 ClpCholeskyDense dense; 2738 // just borrow space 2739 int nDense = numberRows_-firstDense_; 2740 dense.reserveSpace(this,nDense); 2741 dense.solveLong(workDouble_+firstDense_); 2742 for (i=numberRows_-1;i>=firstDense_;i--) { 2743 longDouble value=workDouble_[i]; 2744 int iRow = permute_[i]; 2745 region[iRow]=value; 2746 } 2747 } 2748 for (i=firstDense_-1;i>=0;i--) { 2749 CoinBigIndex offset = indexStart_[i]-choleskyStart_[i]; 2750 longDouble value=workDouble_[i]*diagonal_[i]; 2751 for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) { 2752 int iRow = choleskyRow_[j+offset]; 2753 value -= sparseFactor_[j]*workDouble_[iRow]; 2754 } 2755 workDouble_[i]=value; 2756 int iRow = permute_[i]; 2757 region[iRow]=value; 2758 } 2759 break; 2743 } 2744 for (i=firstDense_-1;i>=0;i--) { 2745 CoinBigIndex offset = indexStart_[i]-choleskyStart_[i]; 2746 CoinWorkDouble value=work[i]*diagonal_[i]; 2747 for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) { 2748 int iRow = choleskyRow_[j+offset]; 2749 value -= sparseFactor_[j]*work[iRow]; 2750 } 2751 work[i]=value; 2752 int iRow = permute_[i]; 2753 region[iRow]=value; 2760 2754 } 2761 2755 } 2762 2756 #endif -
trunk/Clp/src/ClpCholeskyBase.hpp
r1055 r1367 6 6 #include "CoinPragma.hpp" 7 7 #include "CoinFinite.hpp" 8 //#define CLP_LONG_CHOLESKY 0 9 #ifndef CLP_LONG_CHOLESKY 8 10 #define CLP_LONG_CHOLESKY 0 11 #endif 12 /* valid combinations are 13 CLP_LONG_CHOLESKY 0 and COIN_LONG_WORK 0 14 CLP_LONG_CHOLESKY 1 and COIN_LONG_WORK 1 15 CLP_LONG_CHOLESKY 2 and COIN_LONG_WORK 1 16 */ 17 #if COIN_LONG_WORK==0 18 #if CLP_LONG_CHOLESKY>0 19 #define CHOLESKY_BAD_COMBINATION 20 #endif 21 #else 22 #if CLP_LONG_CHOLESKY==0 23 #define CHOLESKY_BAD_COMBINATION 24 #endif 25 #endif 26 #ifdef CHOLESKY_BAD_COMBINATION 27 # warning("Bad combination of CLP_LONG_CHOLESKY and COIN_BIG_DOUBLE/COIN_LONG_WORK"); 28 "Bad combination of CLP_LONG_CHOLESKY and COIN_LONG_WORK" 29 #endif 9 30 #if CLP_LONG_CHOLESKY>1 10 31 typedef long double longDouble; 11 typedef long double longWork; 32 #define CHOL_SMALL_VALUE 1.0e-15 12 33 #elif CLP_LONG_CHOLESKY==1 13 34 typedef double longDouble; 14 typedef long double longWork; 35 #define CHOL_SMALL_VALUE 1.0e-11 15 36 #else 16 37 typedef double longDouble; 17 typedef double longWork; 38 #define CHOL_SMALL_VALUE 1.0e-11 18 39 #endif 19 40 class ClpInterior; … … 46 67 /** Factorize - filling in rowsDropped and returning number dropped. 47 68 If return code negative then out of memory */ 48 virtual int factorize(const double * diagonal, int * rowsDropped) ;69 virtual int factorize(const CoinWorkDouble * diagonal, int * rowsDropped) ; 49 70 /** Uses factorization to solve. */ 50 71 virtual void solve (double * region) ; 51 72 /** Uses factorization to solve. - given as if KKT. 52 73 region1 is rows+columns, region2 is rows */ 53 virtual void solveKKT (double * region1, double * region2, const double * diagonal, 54 double diagonalScaleFactor); 74 virtual void solveKKT (CoinWorkDouble * region1, CoinWorkDouble * region2, const CoinWorkDouble * diagonal, 75 CoinWorkDouble diagonalScaleFactor); 76 #if CLP_LONG_CHOLESKY 77 /** Uses factorization to solve. */ 78 virtual void solve (CoinWorkDouble * region) ; 79 #endif 55 80 //@} 56 81 … … 163 188 */ 164 189 void solve(double * region, int type); 165 void solveLong(longDouble * region, int type);166 190 /// Forms ADAT - returns nonzero if not enough memory 167 191 int preOrder(bool lowerTriangular, bool includeDiagonal, bool doKKT); -
trunk/Clp/src/ClpCholeskyDense.cpp
r1321 r1367 6 6 #include "CoinPragma.hpp" 7 7 #include "CoinHelperFunctions.hpp" 8 #include "ClpHelperFunctions.hpp" 8 9 9 10 #include "ClpInterior.hpp" … … 165 166 CoinZeroN(sparseFactor_,sizeFactor_); 166 167 //perturbation 167 longDouble perturbation=model_->diagonalPerturbation()*model_->diagonalNorm();168 CoinWorkDouble perturbation=model_->diagonalPerturbation()*model_->diagonalNorm(); 168 169 perturbation=perturbation*perturbation; 169 170 if (perturbation>1.0) { … … 172 173 std::cout <<"large perturbation "<<perturbation<<std::endl; 173 174 #endif 174 perturbation= sqrt(perturbation);;175 perturbation=CoinSqrt(perturbation);; 175 176 perturbation=1.0; 176 177 } 177 178 int iRow; 178 179 int newDropped=0; 179 double largest=1.0;180 double smallest=COIN_DBL_MAX;181 longDouble delta2 = model_->delta(); // add delta*delta to diagonal180 CoinWorkDouble largest=1.0; 181 CoinWorkDouble smallest=COIN_DBL_MAX; 182 CoinWorkDouble delta2 = model_->delta(); // add delta*delta to diagonal 182 183 delta2 *= delta2; 183 184 if (!doKKT_) { … … 187 188 const double * diagonalSlack = diagonal + numberColumns; 188 189 // largest in initial matrix 189 double largest2=1.0e-20;190 CoinWorkDouble largest2=1.0e-20; 190 191 for (iRow=0;iRow<numberRows_;iRow++) { 191 192 if (!rowsDropped_[iRow]) { 192 193 CoinBigIndex startRow=rowStart[iRow]; 193 194 CoinBigIndex endRow=rowStart[iRow]+rowLength[iRow]; 194 longDouble diagonalValue = diagonalSlack[iRow]+delta2;195 CoinWorkDouble diagonalValue = diagonalSlack[iRow]+delta2; 195 196 for (CoinBigIndex k=startRow;k<endRow;k++) { 196 197 int iColumn=column[k]; 197 198 CoinBigIndex start=columnStart[iColumn]; 198 199 CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn]; 199 longDouble multiplier = diagonal[iColumn]*elementByRow[k];200 CoinWorkDouble multiplier = diagonal[iColumn]*elementByRow[k]; 200 201 for (CoinBigIndex j=start;j<end;j++) { 201 202 int jRow=row[j]; 202 203 if (!rowsDropped_[jRow]) { 203 204 if (jRow>iRow) { 204 longDouble value=element[j]*multiplier;205 CoinWorkDouble value=element[j]*multiplier; 205 206 work[jRow] += value; 206 207 } else if (jRow==iRow) { 207 longDouble value=element[j]*multiplier;208 CoinWorkDouble value=element[j]*multiplier; 208 209 diagonalValue += value; 209 210 } … … 212 213 } 213 214 for (int j=iRow+1;j<numberRows_;j++) 214 largest2 = CoinMax(largest2, fabs(work[j]));215 largest2 = CoinMax(largest2,CoinAbs(work[j])); 215 216 diagonal_[iRow]=diagonalValue; 216 largest2 = CoinMax(largest2, fabs(diagonalValue));217 largest2 = CoinMax(largest2,CoinAbs(diagonalValue)); 217 218 } else { 218 219 // dropped … … 224 225 //check sizes 225 226 largest2*=1.0e-20; 226 largest = CoinMin(largest2, 1.0e-11);227 largest = CoinMin(largest2,CHOL_SMALL_VALUE); 227 228 int numberDroppedBefore=0; 228 229 for (iRow=0;iRow<numberRows_;iRow++) { … … 231 232 rowsDropped[iRow]=dropped; 232 233 if (!dropped) { 233 longDouble diagonal = diagonal_[iRow];234 CoinWorkDouble diagonal = diagonal_[iRow]; 234 235 if (diagonal>largest2) { 235 236 diagonal_[iRow]=diagonal+perturbation; … … 289 290 if (!quadratic) { 290 291 for (iColumn=0;iColumn<numberColumns;iColumn++) { 291 longDouble value = diagonal[iColumn];292 if ( fabs(value)>1.0e-100) {292 CoinWorkDouble value = diagonal[iColumn]; 293 if (CoinAbs(value)>1.0e-100) { 293 294 value = 1.0/value; 294 largest = CoinMax(largest, fabs(value));295 largest = CoinMax(largest,CoinAbs(value)); 295 296 diagonal_[iColumn] = -value; 296 297 CoinBigIndex start=columnStart[iColumn]; … … 300 301 //sparseFactor_[numberElements++]=element[j]; 301 302 work[row[j]+numberTotal]=element[j]; 302 largest = CoinMax(largest, fabs(element[j]));303 largest = CoinMax(largest,CoinAbs(element[j])); 303 304 } 304 305 } else { … … 315 316 const double * quadraticElement = quadratic->getElements(); 316 317 for (iColumn=0;iColumn<numberColumns;iColumn++) { 317 longDouble value = diagonal[iColumn];318 CoinWorkDouble value = diagonal[iColumn]; 318 319 CoinBigIndex j; 319 if ( fabs(value)>1.0e-100) {320 if (CoinAbs(value)>1.0e-100) { 320 321 value = 1.0/value; 321 322 for (j=columnQuadraticStart[iColumn]; … … 328 329 } 329 330 } 330 largest = CoinMax(largest, fabs(value));331 largest = CoinMax(largest,CoinAbs(value)); 331 332 diagonal_[iColumn] = -value; 332 333 CoinBigIndex start=columnStart[iColumn]; … … 334 335 for (j=start;j<end;j++) { 335 336 work[row[j]+numberTotal]=element[j]; 336 largest = CoinMax(largest, fabs(element[j]));337 largest = CoinMax(largest,CoinAbs(element[j])); 337 338 } 338 339 } else { … … 346 347 // slacks 347 348 for (iColumn=numberColumns;iColumn<numberTotal;iColumn++) { 348 longDouble value = diagonal[iColumn];349 if ( fabs(value)>1.0e-100) {349 CoinWorkDouble value = diagonal[iColumn]; 350 if (CoinAbs(value)>1.0e-100) { 350 351 value = 1.0/value; 351 largest = CoinMax(largest, fabs(value));352 largest = CoinMax(largest,CoinAbs(value)); 352 353 } else { 353 354 value = 1.0e100; … … 364 365 //check sizes 365 366 largest*=1.0e-20; 366 largest = CoinMin(largest, 1.0e-11);367 largest = CoinMin(largest,CHOL_SMALL_VALUE); 367 368 doubleParameters_[10]=CoinMax(1.0e-20,largest); 368 369 integerParameters_[20]=0; … … 387 388 // Should save adjustments in ..R_ 388 389 int n1=0,n2=0; 389 double * primalR = model_->primalR();390 double * dualR = model_->dualR();390 CoinWorkDouble * primalR = model_->primalR(); 391 CoinWorkDouble * dualR = model_->dualR(); 391 392 for (iRow=0;iRow<numberTotal;iRow++) { 392 393 if (rowsDropped2[iRow]) { … … 429 430 CoinMemcpyN(yy,numberRows_,diagonal_); 430 431 int numberDropped=0; 431 double largest=0.0;432 double smallest=COIN_DBL_MAX;432 CoinWorkDouble largest=0.0; 433 CoinWorkDouble smallest=COIN_DBL_MAX; 433 434 double dropValue = doubleParameters_[10]; 434 435 int firstPositive=integerParameters_[34]; … … 441 442 int addOffsetNow = numberRows_-1;; 442 443 longDouble * workNow=sparseFactor_-1+iColumn; 443 double diagonalValue = diagonal_[iColumn];444 CoinWorkDouble diagonalValue = diagonal_[iColumn]; 444 445 for (iRow=0;iRow<iColumn;iRow++) { 445 446 double aj = *workNow; 446 447 addOffsetNow--; 447 448 workNow += addOffsetNow; 448 double multiplier = workDouble_[iRow]; 449 diagonalValue -= aj*aj*multiplier; 449 diagonalValue -= aj*aj*workDouble_[iRow]; 450 450 } 451 451 bool dropColumn=false; … … 776 776 longDouble * diagonal, longDouble * work, int * rowsDropped) 777 777 { 778 longDouble largest=doubleParameters_[3];779 longDouble smallest=doubleParameters_[4];778 CoinWorkDouble largest=doubleParameters_[3]; 779 CoinWorkDouble smallest=doubleParameters_[4]; 780 780 double dropValue = doubleParameters_[10]; 781 781 int firstPositive=integerParameters_[34]; … … 783 783 int numberDropped=0; 784 784 int i, j, k; 785 longDouble t00,temp1;785 CoinWorkDouble t00,temp1; 786 786 longDouble * aa; 787 787 aa = a-BLOCK; … … 790 790 t00 = aa[j]; 791 791 for (k = 0; k < j; ++k) { 792 longDouble multiplier = work[k];792 CoinWorkDouble multiplier = work[k]; 793 793 t00 -= a[j + k * BLOCK] * a[j + k * BLOCK] * multiplier; 794 794 } 795 795 bool dropColumn=false; 796 longDouble useT00=t00;796 CoinWorkDouble useT00=t00; 797 797 if (j+rowOffset<firstPositive) { 798 798 // must be negative … … 831 831 t00=aa[i]; 832 832 for (k = 0; k < j; ++k) { 833 longDouble multiplier = work[k];833 CoinWorkDouble multiplier = work[k]; 834 834 t00 -= a[i + k * BLOCK] * a[j + k * BLOCK] * multiplier; 835 835 } … … 874 874 for (j = 0; j < BLOCK; j +=2) { 875 875 aa+=2*BLOCK; 876 longDouble temp0 = diagonal[j];877 longDouble temp1 = diagonal[j+1];876 CoinWorkDouble temp0 = diagonal[j]; 877 CoinWorkDouble temp1 = diagonal[j+1]; 878 878 for (int i=0;i<BLOCK;i+=2) { 879 longDouble t00=aUnder[i+j*BLOCK];880 longDouble t10=aUnder[i+ BLOCK + j*BLOCK];881 longDouble t01=aUnder[i+1+j*BLOCK];882 longDouble t11=aUnder[i+1+ BLOCK + j*BLOCK];879 CoinWorkDouble t00=aUnder[i+j*BLOCK]; 880 CoinWorkDouble t10=aUnder[i+ BLOCK + j*BLOCK]; 881 CoinWorkDouble t01=aUnder[i+1+j*BLOCK]; 882 CoinWorkDouble t11=aUnder[i+1+ BLOCK + j*BLOCK]; 883 883 for (int k = 0; k < j; ++k) { 884 longDouble multiplier=work[k];885 longDouble au0 = aUnder[i + k * BLOCK]*multiplier;886 longDouble au1 = aUnder[i + 1 + k * BLOCK]*multiplier;887 longDouble at0 = aTri[j + k * BLOCK];888 longDouble at1 = aTri[j + 1 + k * BLOCK];884 CoinWorkDouble multiplier=work[k]; 885 CoinWorkDouble au0 = aUnder[i + k * BLOCK]*multiplier; 886 CoinWorkDouble au1 = aUnder[i + 1 + k * BLOCK]*multiplier; 887 CoinWorkDouble at0 = aTri[j + k * BLOCK]; 888 CoinWorkDouble at1 = aTri[j + 1 + k * BLOCK]; 889 889 t00 -= au0 * at0; 890 890 t10 -= au0 * at1; … … 893 893 } 894 894 t00 *= temp0; 895 longDouble at1 = aTri[j + 1 + j*BLOCK]*work[j];895 CoinWorkDouble at1 = aTri[j + 1 + j*BLOCK]*work[j]; 896 896 t10 -= t00 * at1; 897 897 t01 *= temp0; … … 908 908 for (j = 0; j < BLOCK; j ++) { 909 909 aa+=BLOCK; 910 longDouble temp1 = diagonal[j];910 CoinWorkDouble temp1 = diagonal[j]; 911 911 for (int i=0;i<nUnder;i++) { 912 longDouble t00=aUnder[i+j*BLOCK];912 CoinWorkDouble t00=aUnder[i+j*BLOCK]; 913 913 for (int k = 0; k < j; ++k) { 914 longDouble multiplier=work[k];914 CoinWorkDouble multiplier=work[k]; 915 915 t00 -= aUnder[i + k * BLOCK] * aTri[j + k * BLOCK] * multiplier; 916 916 } … … 937 937 #endif 938 938 int i, j, k; 939 longDouble t00;939 CoinWorkDouble t00; 940 940 longDouble * aa; 941 941 #ifdef BLOCKUNROLL … … 944 944 aa = aTri-2*BLOCK; 945 945 for (j = 0; j < BLOCK; j +=2) { 946 longDouble t00,t01,t10,t11;946 CoinWorkDouble t00,t01,t10,t11; 947 947 aa+=2*BLOCK; 948 948 aUnder2+=2; … … 951 951 t10=aa[j+1+BLOCK]; 952 952 for (k = 0; k < BLOCK; ++k) { 953 longDouble multiplier = work[k];954 longDouble a0 = aUnder2[k * BLOCK];955 longDouble a1 = aUnder2[1 + k * BLOCK];956 longDouble x0 = a0 * multiplier;957 longDouble x1 = a1 * multiplier;953 CoinWorkDouble multiplier = work[k]; 954 CoinWorkDouble a0 = aUnder2[k * BLOCK]; 955 CoinWorkDouble a1 = aUnder2[1 + k * BLOCK]; 956 CoinWorkDouble x0 = a0 * multiplier; 957 CoinWorkDouble x1 = a1 * multiplier; 958 958 t00 -= a0 * x0; 959 959 t01 -= a1 * x0; … … 969 969 t11=aa[i+1+BLOCK]; 970 970 for (k = 0; k < BLOCK; ++k) { 971 longDouble multiplier = work[k];972 longDouble a0 = aUnder2[k * BLOCK]*multiplier;973 longDouble a1 = aUnder2[1 + k * BLOCK]*multiplier;971 CoinWorkDouble multiplier = work[k]; 972 CoinWorkDouble a0 = aUnder2[k * BLOCK]*multiplier; 973 CoinWorkDouble a1 = aUnder2[1 + k * BLOCK]*multiplier; 974 974 t00 -= aUnder[i + k * BLOCK] * a0; 975 975 t01 -= aUnder[i + k * BLOCK] * a1; … … 991 991 t00=aa[i]; 992 992 for (k = 0; k < BLOCK; ++k) { 993 longDouble multiplier = work[k];993 CoinWorkDouble multiplier = work[k]; 994 994 t00 -= aUnder[i + k * BLOCK] * aUnder[j + k * BLOCK]*multiplier; 995 995 } … … 1034 1034 aa+=2*BLOCK; 1035 1035 for (i=0;i< BLOCK;i+=2) { 1036 longDouble t00=aa[i+0*BLOCK];1037 longDouble t10=aa[i+1*BLOCK];1038 longDouble t01=aa[i+1+0*BLOCK];1039 longDouble t11=aa[i+1+1*BLOCK];1036 CoinWorkDouble t00=aa[i+0*BLOCK]; 1037 CoinWorkDouble t10=aa[i+1*BLOCK]; 1038 CoinWorkDouble t01=aa[i+1+0*BLOCK]; 1039 CoinWorkDouble t11=aa[i+1+1*BLOCK]; 1040 1040 for (k=0;k<BLOCK;k++) { 1041 longDouble multiplier = work[k];1042 longDouble a00=aUnder[i+k*BLOCK]*multiplier;1043 longDouble a01=aUnder[i+1+k*BLOCK]*multiplier;1041 CoinWorkDouble multiplier = work[k]; 1042 CoinWorkDouble a00=aUnder[i+k*BLOCK]*multiplier; 1043 CoinWorkDouble a01=aUnder[i+1+k*BLOCK]*multiplier; 1044 1044 t00 -= a00 * above[j + 0 + k * BLOCK]; 1045 1045 t10 -= a00 * above[j + 1 + k * BLOCK]; … … 1057 1057 aa+=4*BLOCK; 1058 1058 for (i=0;i< BLOCK;i+=4) { 1059 longDouble t00=aa[i+0+0*BLOCK];1060 longDouble t10=aa[i+0+1*BLOCK];1061 longDouble t20=aa[i+0+2*BLOCK];1062 longDouble t30=aa[i+0+3*BLOCK];1063 longDouble t01=aa[i+1+0*BLOCK];1064 longDouble t11=aa[i+1+1*BLOCK];1065 longDouble t21=aa[i+1+2*BLOCK];1066 longDouble t31=aa[i+1+3*BLOCK];1067 longDouble t02=aa[i+2+0*BLOCK];1068 longDouble t12=aa[i+2+1*BLOCK];1069 longDouble t22=aa[i+2+2*BLOCK];1070 longDouble t32=aa[i+2+3*BLOCK];1071 longDouble t03=aa[i+3+0*BLOCK];1072 longDouble t13=aa[i+3+1*BLOCK];1073 longDouble t23=aa[i+3+2*BLOCK];1074 longDouble t33=aa[i+3+3*BLOCK];1059 CoinWorkDouble t00=aa[i+0+0*BLOCK]; 1060 CoinWorkDouble t10=aa[i+0+1*BLOCK]; 1061 CoinWorkDouble t20=aa[i+0+2*BLOCK]; 1062 CoinWorkDouble t30=aa[i+0+3*BLOCK]; 1063 CoinWorkDouble t01=aa[i+1+0*BLOCK]; 1064 CoinWorkDouble t11=aa[i+1+1*BLOCK]; 1065 CoinWorkDouble t21=aa[i+1+2*BLOCK]; 1066 CoinWorkDouble t31=aa[i+1+3*BLOCK]; 1067 CoinWorkDouble t02=aa[i+2+0*BLOCK]; 1068 CoinWorkDouble t12=aa[i+2+1*BLOCK]; 1069 CoinWorkDouble t22=aa[i+2+2*BLOCK]; 1070 CoinWorkDouble t32=aa[i+2+3*BLOCK]; 1071 CoinWorkDouble t03=aa[i+3+0*BLOCK]; 1072 CoinWorkDouble t13=aa[i+3+1*BLOCK]; 1073 CoinWorkDouble t23=aa[i+3+2*BLOCK]; 1074 CoinWorkDouble t33=aa[i+3+3*BLOCK]; 1075 1075 const longDouble * COIN_RESTRICT aUnderNow = aUnder+i; 1076 1076 const longDouble * COIN_RESTRICT aboveNow = above+j; 1077 1077 for (k=0;k<BLOCK;k++) { 1078 longDouble multiplier = work[k];1079 longDouble a00=aUnderNow[0]*multiplier;1080 longDouble a01=aUnderNow[1]*multiplier;1081 longDouble a02=aUnderNow[2]*multiplier;1082 longDouble a03=aUnderNow[3]*multiplier;1078 CoinWorkDouble multiplier = work[k]; 1079 CoinWorkDouble a00=aUnderNow[0]*multiplier; 1080 CoinWorkDouble a01=aUnderNow[1]*multiplier; 1081 CoinWorkDouble a02=aUnderNow[2]*multiplier; 1082 CoinWorkDouble a03=aUnderNow[3]*multiplier; 1083 1083 t00 -= a00 * aboveNow[0]; 1084 1084 t10 -= a00 * aboveNow[1]; … … 1125 1125 aa+=4*BLOCK; 1126 1126 for (i=0;i< n;i+=2) { 1127 longDouble t00=aa[i+0*BLOCK];1128 longDouble t10=aa[i+1*BLOCK];1129 longDouble t20=aa[i+2*BLOCK];1130 longDouble t30=aa[i+3*BLOCK];1131 longDouble t01=aa[i+1+0*BLOCK];1132 longDouble t11=aa[i+1+1*BLOCK];1133 longDouble t21=aa[i+1+2*BLOCK];1134 longDouble t31=aa[i+1+3*BLOCK];1127 CoinWorkDouble t00=aa[i+0*BLOCK]; 1128 CoinWorkDouble t10=aa[i+1*BLOCK]; 1129 CoinWorkDouble t20=aa[i+2*BLOCK]; 1130 CoinWorkDouble t30=aa[i+3*BLOCK]; 1131 CoinWorkDouble t01=aa[i+1+0*BLOCK]; 1132 CoinWorkDouble t11=aa[i+1+1*BLOCK]; 1133 CoinWorkDouble t21=aa[i+1+2*BLOCK]; 1134 CoinWorkDouble t31=aa[i+1+3*BLOCK]; 1135 1135 const longDouble * COIN_RESTRICT aUnderNow = aUnder+i; 1136 1136 const longDouble * COIN_RESTRICT aboveNow = above+j; 1137 1137 for (k=0;k<BLOCK;k++) { 1138 longDouble multiplier = work[k];1139 longDouble a00=aUnderNow[0]*multiplier;1140 longDouble a01=aUnderNow[1]*multiplier;1138 CoinWorkDouble multiplier = work[k]; 1139 CoinWorkDouble a00=aUnderNow[0]*multiplier; 1140 CoinWorkDouble a01=aUnderNow[1]*multiplier; 1141 1141 t00 -= a00 * aboveNow[0]; 1142 1142 t10 -= a00 * aboveNow[1]; … … 1160 1160 } 1161 1161 if (odd) { 1162 longDouble t0=aa[n+0*BLOCK];1163 longDouble t1=aa[n+1*BLOCK];1164 longDouble t2=aa[n+2*BLOCK];1165 longDouble t3=aa[n+3*BLOCK];1166 longDouble a0;1162 CoinWorkDouble t0=aa[n+0*BLOCK]; 1163 CoinWorkDouble t1=aa[n+1*BLOCK]; 1164 CoinWorkDouble t2=aa[n+2*BLOCK]; 1165 CoinWorkDouble t3=aa[n+3*BLOCK]; 1166 CoinWorkDouble a0; 1167 1167 for (k=0;k<BLOCK;k++) { 1168 1168 a0=aUnder[n+k*BLOCK]*work[k]; … … 1184 1184 aa+=BLOCK; 1185 1185 for (i=0;i< nUnder;i++) { 1186 longDouble t00=aa[i+0*BLOCK];1186 CoinWorkDouble t00=aa[i+0*BLOCK]; 1187 1187 for (k=0;k<BLOCK;k++) { 1188 longDouble a00=aUnder[i+k*BLOCK]*work[k];1188 CoinWorkDouble a00=aUnder[i+k*BLOCK]*work[k]; 1189 1189 t00 -= a00 * above[j + k * BLOCK]; 1190 1190 } … … 1293 1293 if (numberRows_<200) { 1294 1294 for (int i=0;i<numberRows_;i++) { 1295 assert( fabs(region[i]-region2[i])<1.0e-3);1295 assert(CoinAbs(region[i]-region2[i])<1.0e-3); 1296 1296 } 1297 1297 delete [] region2; … … 1304 1304 { 1305 1305 int j, k; 1306 longDouble t00;1306 CoinWorkDouble t00; 1307 1307 for (j = 0; j < n; j ++) { 1308 1308 t00 = region[j]; … … 1322 1322 if (n==BLOCK) { 1323 1323 for (k = 0; k < BLOCK; k+=4) { 1324 longDouble t0 = region2[0];1325 longDouble t1 = region2[1];1326 longDouble t2 = region2[2];1327 longDouble t3 = region2[3];1324 CoinWorkDouble t0 = region2[0]; 1325 CoinWorkDouble t1 = region2[1]; 1326 CoinWorkDouble t2 = region2[2]; 1327 CoinWorkDouble t3 = region2[3]; 1328 1328 t0 -= region[0] * a[0 + 0 * BLOCK]; 1329 1329 t1 -= region[0] * a[1 + 0 * BLOCK]; … … 1416 1416 #endif 1417 1417 for (k = 0; k < n; ++k) { 1418 longDouble t00 = region2[k];1418 CoinWorkDouble t00 = region2[k]; 1419 1419 for (j = 0; j < BLOCK; j ++) { 1420 1420 t00 -= region[j] * a[k + j * BLOCK]; … … 1431 1431 { 1432 1432 int j, k; 1433 longDouble t00;1433 CoinWorkDouble t00; 1434 1434 for (j = n-1; j >=0; j --) { 1435 1435 t00 = region[j]; … … 1449 1449 if (n==BLOCK) { 1450 1450 for (j = 0; j < BLOCK; j +=4) { 1451 longDouble t0 = region[0];1452 longDouble t1 = region[1];1453 longDouble t2 = region[2];1454 longDouble t3 = region[3];1451 CoinWorkDouble t0 = region[0]; 1452 CoinWorkDouble t1 = region[1]; 1453 CoinWorkDouble t2 = region[2]; 1454 CoinWorkDouble t3 = region[3]; 1455 1455 t0 -= region2[0] * a[0 + 0*BLOCK]; 1456 1456 t1 -= region2[0] * a[0 + 1*BLOCK]; … … 1544 1544 #endif 1545 1545 for (j = 0; j < BLOCK; j ++) { 1546 longDouble t00 = region[j];1546 CoinWorkDouble t00 = region[j]; 1547 1547 for (k = 0; k < n; ++k) { 1548 1548 t00 -= region2[k] * a[k + j * BLOCK]; … … 1556 1556 /* Uses factorization to solve. */ 1557 1557 void 1558 ClpCholeskyDense::solveLong ( longDouble * region)1558 ClpCholeskyDense::solveLong (CoinWorkDouble * region) 1559 1559 { 1560 1560 int numberBlocks = (numberRows_+BLOCK-1)>>BLOCKSHIFT; … … 1619 1619 // Forward part of solve 1 1620 1620 void 1621 ClpCholeskyDense::solveF1Long(longDouble * a,int n, longDouble * region)1621 ClpCholeskyDense::solveF1Long(longDouble * a,int n,CoinWorkDouble * region) 1622 1622 { 1623 1623 int j, k; 1624 longDouble t00;1624 CoinWorkDouble t00; 1625 1625 for (j = 0; j < n; j ++) { 1626 1626 t00 = region[j]; … … 1634 1634 // Forward part of solve 2 1635 1635 void 1636 ClpCholeskyDense::solveF2Long(longDouble * a,int n, longDouble * region, longDouble * region2)1636 ClpCholeskyDense::solveF2Long(longDouble * a,int n,CoinWorkDouble * region, CoinWorkDouble * region2) 1637 1637 { 1638 1638 int j, k; … … 1640 1640 if (n==BLOCK) { 1641 1641 for (k = 0; k < BLOCK; k+=4) { 1642 longDouble t0 = region2[0];1643 longDouble t1 = region2[1];1644 longDouble t2 = region2[2];1645 longDouble t3 = region2[3];1642 CoinWorkDouble t0 = region2[0]; 1643 CoinWorkDouble t1 = region2[1]; 1644 CoinWorkDouble t2 = region2[2]; 1645 CoinWorkDouble t3 = region2[3]; 1646 1646 t0 -= region[0] * a[0 + 0 * BLOCK]; 1647 1647 t1 -= region[0] * a[1 + 0 * BLOCK]; … … 1734 1734 #endif 1735 1735 for (k = 0; k < n; ++k) { 1736 longDouble t00 = region2[k];1736 CoinWorkDouble t00 = region2[k]; 1737 1737 for (j = 0; j < BLOCK; j ++) { 1738 1738 t00 -= region[j] * a[k + j * BLOCK]; … … 1746 1746 // Backward part of solve 1 1747 1747 void 1748 ClpCholeskyDense::solveB1Long(longDouble * a,int n, longDouble * region)1748 ClpCholeskyDense::solveB1Long(longDouble * a,int n,CoinWorkDouble * region) 1749 1749 { 1750 1750 int j, k; 1751 longDouble t00;1751 CoinWorkDouble t00; 1752 1752 for (j = n-1; j >=0; j --) { 1753 1753 t00 = region[j]; … … 1761 1761 // Backward part of solve 2 1762 1762 void 1763 ClpCholeskyDense::solveB2Long(longDouble * a,int n, longDouble * region, longDouble * region2)1763 ClpCholeskyDense::solveB2Long(longDouble * a,int n,CoinWorkDouble * region, CoinWorkDouble * region2) 1764 1764 { 1765 1765 int j, k; … … 1767 1767 if (n==BLOCK) { 1768 1768 for (j = 0; j < BLOCK; j +=4) { 1769 longDouble t0 = region[0];1770 longDouble t1 = region[1];1771 longDouble t2 = region[2];1772 longDouble t3 = region[3];1769 CoinWorkDouble t0 = region[0]; 1770 CoinWorkDouble t1 = region[1]; 1771 CoinWorkDouble t2 = region[2]; 1772 CoinWorkDouble t3 = region[3]; 1773 1773 t0 -= region2[0] * a[0 + 0*BLOCK]; 1774 1774 t1 -= region2[0] * a[0 + 1*BLOCK]; … … 1862 1862 #endif 1863 1863 for (j = 0; j < BLOCK; j ++) { 1864 longDouble t00 = region[j];1864 CoinWorkDouble t00 = region[j]; 1865 1865 for (k = 0; k < n; ++k) { 1866 1866 t00 -= region2[k] * a[k + j * BLOCK]; … … 1874 1874 /* Uses factorization to solve. */ 1875 1875 void 1876 ClpCholeskyDense::solveLongWork ( longWork* region)1876 ClpCholeskyDense::solveLongWork (CoinWorkDouble * region) 1877 1877 { 1878 1878 int numberBlocks = (numberRows_+BLOCK-1)>>BLOCKSHIFT; … … 1937 1937 // Forward part of solve 1 1938 1938 void 1939 ClpCholeskyDense::solveF1LongWork(longDouble * a,int n, longWork* region)1939 ClpCholeskyDense::solveF1LongWork(longDouble * a,int n,CoinWorkDouble * region) 1940 1940 { 1941 1941 int j, k; 1942 longWorkt00;1942 CoinWorkDouble t00; 1943 1943 for (j = 0; j < n; j ++) { 1944 1944 t00 = region[j]; … … 1952 1952 // Forward part of solve 2 1953 1953 void 1954 ClpCholeskyDense::solveF2LongWork(longDouble * a,int n, longWork * region, longWork* region2)1954 ClpCholeskyDense::solveF2LongWork(longDouble * a,int n,CoinWorkDouble * region, CoinWorkDouble * region2) 1955 1955 { 1956 1956 int j, k; … … 1958 1958 if (n==BLOCK) { 1959 1959 for (k = 0; k < BLOCK; k+=4) { 1960 longWorkt0 = region2[0];1961 longWorkt1 = region2[1];1962 longWorkt2 = region2[2];1963 longWorkt3 = region2[3];1960 CoinWorkDouble t0 = region2[0]; 1961 CoinWorkDouble t1 = region2[1]; 1962 CoinWorkDouble t2 = region2[2]; 1963 CoinWorkDouble t3 = region2[3]; 1964 1964 t0 -= region[0] * a[0 + 0 * BLOCK]; 1965 1965 t1 -= region[0] * a[1 + 0 * BLOCK]; … … 2052 2052 #endif 2053 2053 for (k = 0; k < n; ++k) { 2054 longWorkt00 = region2[k];2054 CoinWorkDouble t00 = region2[k]; 2055 2055 for (j = 0; j < BLOCK; j ++) { 2056 2056 t00 -= region[j] * a[k + j * BLOCK]; … … 2064 2064 // Backward part of solve 1 2065 2065 void 2066 ClpCholeskyDense::solveB1LongWork(longDouble * a,int n, longWork* region)2066 ClpCholeskyDense::solveB1LongWork(longDouble * a,int n,CoinWorkDouble * region) 2067 2067 { 2068 2068 int j, k; 2069 longWorkt00;2069 CoinWorkDouble t00; 2070 2070 for (j = n-1; j >=0; j --) { 2071 2071 t00 = region[j]; … … 2079 2079 // Backward part of solve 2 2080 2080 void 2081 ClpCholeskyDense::solveB2LongWork(longDouble * a,int n, longWork * region, longWork* region2)2081 ClpCholeskyDense::solveB2LongWork(longDouble * a,int n,CoinWorkDouble * region, CoinWorkDouble * region2) 2082 2082 { 2083 2083 int j, k; … … 2085 2085 if (n==BLOCK) { 2086 2086 for (j = 0; j < BLOCK; j +=4) { 2087 longWorkt0 = region[0];2088 longWorkt1 = region[1];2089 longWorkt2 = region[2];2090 longWorkt3 = region[3];2087 CoinWorkDouble t0 = region[0]; 2088 CoinWorkDouble t1 = region[1]; 2089 CoinWorkDouble t2 = region[2]; 2090 CoinWorkDouble t3 = region[3]; 2091 2091 t0 -= region2[0] * a[0 + 0*BLOCK]; 2092 2092 t1 -= region2[0] * a[0 + 1*BLOCK]; … … 2180 2180 #endif 2181 2181 for (j = 0; j < BLOCK; j ++) { 2182 longWorkt00 = region[j];2182 CoinWorkDouble t00 = region[j]; 2183 2183 for (k = 0; k < n; ++k) { 2184 2184 t00 -= region2[k] * a[k + j * BLOCK]; -
trunk/Clp/src/ClpCholeskyDense.hpp
r1321 r1367 89 89 void solveB2(longDouble * a,int n,double * region,double * region2); 90 90 /** Uses factorization to solve. */ 91 void solveLong ( longDouble * region) ;91 void solveLong (CoinWorkDouble * region) ; 92 92 /// Forward part of solve 93 void solveF1Long(longDouble * a,int n, longDouble * region);94 void solveF2Long(longDouble * a,int n, longDouble * region,longDouble * region2);93 void solveF1Long(longDouble * a,int n,CoinWorkDouble * region); 94 void solveF2Long(longDouble * a,int n,CoinWorkDouble * region,CoinWorkDouble * region2); 95 95 /// Backward part of solve 96 void solveB1Long(longDouble * a,int n, longDouble * region);97 void solveB2Long(longDouble * a,int n, longDouble * region,longDouble * region2);96 void solveB1Long(longDouble * a,int n,CoinWorkDouble * region); 97 void solveB2Long(longDouble * a,int n,CoinWorkDouble * region,CoinWorkDouble * region2); 98 98 /** Uses factorization to solve. */ 99 void solveLongWork ( longWork* region) ;99 void solveLongWork (CoinWorkDouble * region) ; 100 100 /// Forward part of solve 101 void solveF1LongWork(longDouble * a,int n, longWork* region);102 void solveF2LongWork(longDouble * a,int n, longWork * region,longWork* region2);101 void solveF1LongWork(longDouble * a,int n,CoinWorkDouble * region); 102 void solveF2LongWork(longDouble * a,int n,CoinWorkDouble * region,CoinWorkDouble * region2); 103 103 /// Backward part of solve 104 void solveB1LongWork(longDouble * a,int n, longWork* region);105 void solveB2LongWork(longDouble * a,int n, longWork * region,longWork* region2);104 void solveB1LongWork(longDouble * a,int n,CoinWorkDouble * region); 105 void solveB2LongWork(longDouble * a,int n,CoinWorkDouble * region,CoinWorkDouble * region2); 106 106 int bNumber(const longDouble * array,int &, int&); 107 107 /// A -
trunk/Clp/src/ClpHelperFunctions.cpp
r1185 r1367 112 112 } 113 113 } 114 #if COIN_LONG_WORK 115 // For long double versions 116 CoinWorkDouble 117 maximumAbsElement(const CoinWorkDouble * region, int size) 118 { 119 int i; 120 CoinWorkDouble maxValue=0.0; 121 for (i=0;i<size;i++) 122 maxValue = CoinMax(maxValue,CoinAbs(region[i])); 123 return maxValue; 124 } 125 void 126 setElements(CoinWorkDouble * region, int size, CoinWorkDouble value) 127 { 128 int i; 129 for (i=0;i<size;i++) 130 region[i]=value; 131 } 132 void 133 multiplyAdd(const CoinWorkDouble * region1, int size, CoinWorkDouble multiplier1, 134 CoinWorkDouble * region2, CoinWorkDouble multiplier2) 135 { 136 int i; 137 if (multiplier1==1.0) { 138 if (multiplier2==1.0) { 139 for (i=0;i<size;i++) 140 region2[i] = region1[i] + region2[i]; 141 } else if (multiplier2==-1.0) { 142 for (i=0;i<size;i++) 143 region2[i] = region1[i] - region2[i]; 144 } else if (multiplier2==0.0) { 145 for (i=0;i<size;i++) 146 region2[i] = region1[i] ; 147 } else { 148 for (i=0;i<size;i++) 149 region2[i] = region1[i] + multiplier2*region2[i]; 150 } 151 } else if (multiplier1==-1.0) { 152 if (multiplier2==1.0) { 153 for (i=0;i<size;i++) 154 region2[i] = -region1[i] + region2[i]; 155 } else if (multiplier2==-1.0) { 156 for (i=0;i<size;i++) 157 region2[i] = -region1[i] - region2[i]; 158 } else if (multiplier2==0.0) { 159 for (i=0;i<size;i++) 160 region2[i] = -region1[i] ; 161 } else { 162 for (i=0;i<size;i++) 163 region2[i] = -region1[i] + multiplier2*region2[i]; 164 } 165 } else if (multiplier1==0.0) { 166 if (multiplier2==1.0) { 167 // nothing to do 168 } else if (multiplier2==-1.0) { 169 for (i=0;i<size;i++) 170 region2[i] = -region2[i]; 171 } else if (multiplier2==0.0) { 172 for (i=0;i<size;i++) 173 region2[i] = 0.0; 174 } else { 175 for (i=0;i<size;i++) 176 region2[i] = multiplier2*region2[i]; 177 } 178 } else { 179 if (multiplier2==1.0) { 180 for (i=0;i<size;i++) 181 region2[i] = multiplier1*region1[i] + region2[i]; 182 } else if (multiplier2==-1.0) { 183 for (i=0;i<size;i++) 184 region2[i] = multiplier1*region1[i] - region2[i]; 185 } else if (multiplier2==0.0) { 186 for (i=0;i<size;i++) 187 region2[i] = multiplier1*region1[i] ; 188 } else { 189 for (i=0;i<size;i++) 190 region2[i] = multiplier1*region1[i] + multiplier2*region2[i]; 191 } 192 } 193 } 194 CoinWorkDouble 195 innerProduct(const CoinWorkDouble * region1, int size, const CoinWorkDouble * region2) 196 { 197 int i; 198 CoinWorkDouble value=0.0; 199 for (i=0;i<size;i++) 200 value += region1[i]*region2[i]; 201 return value; 202 } 203 void 204 getNorms(const CoinWorkDouble * region, int size, CoinWorkDouble & norm1, CoinWorkDouble & norm2) 205 { 206 norm1 = 0.0; 207 norm2 = 0.0; 208 int i; 209 for (i=0;i<size;i++) { 210 norm2 += region[i]*region[i]; 211 norm1 = CoinMax(norm1,CoinAbs(region[i])); 212 } 213 } 214 #endif 114 215 #ifdef DEBUG_MEMORY 115 216 #include <malloc.h> -
trunk/Clp/src/ClpHelperFunctions.hpp
r754 r1367 18 18 double innerProduct(const double * region1, int size, const double * region2); 19 19 void getNorms(const double * region, int size, double & norm1, double & norm2); 20 #if COIN_LONG_WORK 21 // For long double versions 22 CoinWorkDouble maximumAbsElement(const CoinWorkDouble * region, int size); 23 void setElements(CoinWorkDouble * region, int size, CoinWorkDouble value); 24 void multiplyAdd(const CoinWorkDouble * region1, int size, CoinWorkDouble multiplier1, 25 CoinWorkDouble * region2, CoinWorkDouble multiplier2); 26 CoinWorkDouble innerProduct(const CoinWorkDouble * region1, int size, const CoinWorkDouble * region2); 27 void getNorms(const CoinWorkDouble * region, int size, CoinWorkDouble & norm1, CoinWorkDouble & norm2); 28 inline void 29 CoinMemcpyN(const double * from, const int size, CoinWorkDouble * to) 30 { 31 for (int i=0;i<size;i++) 32 to[i]=from[i]; 33 } 34 inline void 35 CoinMemcpyN(const CoinWorkDouble * from, const int size, double * to) 36 { 37 for (int i=0;i<size;i++) 38 to[i]=from[i]; 39 } 40 inline CoinWorkDouble 41 CoinMax(const CoinWorkDouble x1, const double x2) 42 { 43 return (x1 > x2) ? x1 : x2; 44 } 45 inline CoinWorkDouble 46 CoinMax(double x1, const CoinWorkDouble x2) 47 { 48 return (x1 > x2) ? x1 : x2; 49 } 50 inline CoinWorkDouble 51 CoinMin(const CoinWorkDouble x1, const double x2) 52 { 53 return (x1 < x2) ? x1 : x2; 54 } 55 inline CoinWorkDouble 56 CoinMin(double x1, const CoinWorkDouble x2) 57 { 58 return (x1 < x2) ? x1 : x2; 59 } 60 inline CoinWorkDouble CoinSqrt(CoinWorkDouble x) 61 { 62 return sqrtl(x); 63 } 64 #else 65 inline double CoinSqrt(double x) 66 { 67 return sqrt(x); 68 } 69 #endif 20 70 21 71 /// Following only included if ClpPdco defined -
trunk/Clp/src/ClpInterior.cpp
r1321 r1367 109 109 algorithm_(-1) 110 110 { 111 memset(historyInfeasibility_,0,LENGTH_HISTORY*sizeof( double));111 memset(historyInfeasibility_,0,LENGTH_HISTORY*sizeof(CoinWorkDouble)); 112 112 solveType_=3; // say interior based life form 113 113 cholesky_ = new ClpCholeskyDense(); // put in placeholder … … 198 198 algorithm_(-1) 199 199 { 200 memset(historyInfeasibility_,0,LENGTH_HISTORY*sizeof( double));200 memset(historyInfeasibility_,0,LENGTH_HISTORY*sizeof(CoinWorkDouble)); 201 201 solveType_=3; // say interior based life form 202 202 cholesky_= new ClpCholeskyDense(); … … 351 351 algorithm_(-1) 352 352 { 353 memset(historyInfeasibility_,0,LENGTH_HISTORY*sizeof( double));353 memset(historyInfeasibility_,0,LENGTH_HISTORY*sizeof(CoinWorkDouble)); 354 354 solveType_=3; // say interior based life form 355 355 cholesky_= new ClpCholeskyDense(); … … 537 537 int nTotal = numberRows_+numberColumns_; 538 538 delete [] solution_; 539 solution_ = new double[nTotal];539 solution_ = new CoinWorkDouble[nTotal]; 540 540 CoinMemcpyN(columnActivity_, numberColumns_,solution_); 541 541 CoinMemcpyN(rowActivity_, numberRows_,solution_+numberColumns_); 542 542 delete [] cost_; 543 cost_ = new double[nTotal];543 cost_ = new CoinWorkDouble[nTotal]; 544 544 int i; 545 double direction = optimizationDirection_*objectiveScale_;545 CoinWorkDouble direction = optimizationDirection_*objectiveScale_; 546 546 // direction is actually scale out not scale in 547 547 if (direction) … … 550 550 for (i=0;i<numberColumns_;i++) 551 551 cost_[i] = direction*obj[i]; 552 memset(cost_+numberColumns_,0,numberRows_*sizeof( double));552 memset(cost_+numberColumns_,0,numberRows_*sizeof(CoinWorkDouble)); 553 553 // do scaling if needed 554 554 if (scalingFlag_>0&&!rowScale_) { … … 558 558 delete [] lower_; 559 559 delete [] upper_; 560 lower_ = new double[nTotal];561 upper_ = new double[nTotal];560 lower_ = new CoinWorkDouble[nTotal]; 561 upper_ = new CoinWorkDouble[nTotal]; 562 562 rowLowerWork_ = lower_+numberColumns_; 563 563 columnLowerWork_ = lower_; … … 587 587 if(rowScale_) { 588 588 for (i=0;i<numberColumns_;i++) { 589 double multiplier = rhsScale_/columnScale_[i];589 CoinWorkDouble multiplier = rhsScale_/columnScale_[i]; 590 590 cost_[i] *= columnScale_[i]; 591 591 if (columnLowerWork_[i]>-1.0e50) … … 596 596 } 597 597 for (i=0;i<numberRows_;i++) { 598 double multiplier = rhsScale_*rowScale_[i];598 CoinWorkDouble multiplier = rhsScale_*rowScale_[i]; 599 599 if (rowLowerWork_[i]>-1.0e50) 600 600 rowLowerWork_[i] *= multiplier; … … 612 612 } 613 613 assert (!errorRegion_); 614 errorRegion_ = new double [numberRows_];614 errorRegion_ = new CoinWorkDouble [numberRows_]; 615 615 assert (!rhsFixRegion_); 616 rhsFixRegion_ = new double [numberRows_];616 rhsFixRegion_ = new CoinWorkDouble [numberRows_]; 617 617 assert (!deltaY_); 618 deltaY_ = new double [numberRows_];618 deltaY_ = new CoinWorkDouble [numberRows_]; 619 619 CoinZeroN(deltaY_,numberRows_); 620 620 assert (!upperSlack_); 621 upperSlack_ = new double [nTotal];621 upperSlack_ = new CoinWorkDouble [nTotal]; 622 622 assert (!lowerSlack_); 623 lowerSlack_ = new double [nTotal];623 lowerSlack_ = new CoinWorkDouble [nTotal]; 624 624 assert (!diagonal_); 625 diagonal_ = new double [nTotal];625 diagonal_ = new CoinWorkDouble [nTotal]; 626 626 assert (!deltaX_); 627 deltaX_ = new double [nTotal];627 deltaX_ = new CoinWorkDouble [nTotal]; 628 628 CoinZeroN(deltaX_,nTotal); 629 629 assert (!deltaZ_); 630 deltaZ_ = new double [nTotal];630 deltaZ_ = new CoinWorkDouble [nTotal]; 631 631 CoinZeroN(deltaZ_,nTotal); 632 632 assert (!deltaW_); 633 deltaW_ = new double [nTotal];633 deltaW_ = new CoinWorkDouble [nTotal]; 634 634 CoinZeroN(deltaW_,nTotal); 635 635 assert (!deltaSU_); 636 deltaSU_ = new double [nTotal];636 deltaSU_ = new CoinWorkDouble [nTotal]; 637 637 CoinZeroN(deltaSU_,nTotal); 638 638 assert (!deltaSL_); 639 deltaSL_ = new double [nTotal];639 deltaSL_ = new CoinWorkDouble [nTotal]; 640 640 CoinZeroN(deltaSL_,nTotal); 641 641 assert (!primalR_); … … 643 643 // create arrays if we are doing KKT 644 644 if (cholesky_->type()>=20) { 645 primalR_ = new double [nTotal];645 primalR_ = new CoinWorkDouble [nTotal]; 646 646 CoinZeroN(primalR_,nTotal); 647 dualR_ = new double [numberRows_];647 dualR_ = new CoinWorkDouble [numberRows_]; 648 648 CoinZeroN(dualR_,numberRows_); 649 649 } 650 650 assert (!rhsB_); 651 rhsB_ = new double [numberRows_];651 rhsB_ = new CoinWorkDouble [numberRows_]; 652 652 CoinZeroN(rhsB_,numberRows_); 653 653 assert (!rhsU_); 654 rhsU_ = new double [nTotal];654 rhsU_ = new CoinWorkDouble [nTotal]; 655 655 CoinZeroN(rhsU_,nTotal); 656 656 assert (!rhsL_); 657 rhsL_ = new double [nTotal];657 rhsL_ = new CoinWorkDouble [nTotal]; 658 658 CoinZeroN(rhsL_,nTotal); 659 659 assert (!rhsZ_); 660 rhsZ_ = new double [nTotal];660 rhsZ_ = new CoinWorkDouble [nTotal]; 661 661 CoinZeroN(rhsZ_,nTotal); 662 662 assert (!rhsW_); 663 rhsW_ = new double [nTotal];663 rhsW_ = new CoinWorkDouble [nTotal]; 664 664 CoinZeroN(rhsW_,nTotal); 665 665 assert (!rhsC_); 666 rhsC_ = new double [nTotal];666 rhsC_ = new CoinWorkDouble [nTotal]; 667 667 CoinZeroN(rhsC_,nTotal); 668 668 assert (!workArray_); 669 workArray_ = new double [nTotal];669 workArray_ = new CoinWorkDouble [nTotal]; 670 670 CoinZeroN(workArray_,nTotal); 671 671 assert (!zVec_); 672 zVec_ = new double [nTotal];672 zVec_ = new CoinWorkDouble [nTotal]; 673 673 CoinZeroN(zVec_,nTotal); 674 674 assert (!wVec_); 675 wVec_ = new double [nTotal];675 wVec_ = new CoinWorkDouble [nTotal]; 676 676 CoinZeroN(wVec_,nTotal); 677 677 assert (!dj_); 678 dj_ = new double [nTotal];678 dj_ = new CoinWorkDouble [nTotal]; 679 679 if (!status_) 680 680 status_ = new unsigned char [numberRows_+numberColumns_]; … … 687 687 int i; 688 688 if (optimizationDirection_!=1.0||objectiveScale_!=1.0) { 689 double scaleC = optimizationDirection_/objectiveScale_;689 CoinWorkDouble scaleC = optimizationDirection_/objectiveScale_; 690 690 // and modify all dual signs 691 691 for (i=0;i<numberColumns_;i++) … … 695 695 } 696 696 if (rowScale_) { 697 double scaleR = 1.0/rhsScale_;697 CoinWorkDouble scaleR = 1.0/rhsScale_; 698 698 for (i=0;i<numberColumns_;i++) { 699 double scaleFactor = columnScale_[i];700 double valueScaled = columnActivity_[i];699 CoinWorkDouble scaleFactor = columnScale_[i]; 700 CoinWorkDouble valueScaled = columnActivity_[i]; 701 701 columnActivity_[i] = valueScaled*scaleFactor*scaleR; 702 double valueScaledDual = reducedCost_[i];702 CoinWorkDouble valueScaledDual = reducedCost_[i]; 703 703 reducedCost_[i] = valueScaledDual/scaleFactor; 704 704 } 705 705 for (i=0;i<numberRows_;i++) { 706 double scaleFactor = rowScale_[i];707 double valueScaled = rowActivity_[i];706 CoinWorkDouble scaleFactor = rowScale_[i]; 707 CoinWorkDouble valueScaled = rowActivity_[i]; 708 708 rowActivity_[i] = (valueScaled*scaleR)/scaleFactor; 709 double valueScaledDual = dual_[i];709 CoinWorkDouble valueScaledDual = dual_[i]; 710 710 dual_[i] = valueScaledDual*scaleFactor; 711 711 } 712 712 } else if (rhsScale_!=1.0) { 713 double scaleR = 1.0/rhsScale_;713 CoinWorkDouble scaleR = 1.0/rhsScale_; 714 714 for (i=0;i<numberColumns_;i++) { 715 double valueScaled = columnActivity_[i];715 CoinWorkDouble valueScaled = columnActivity_[i]; 716 716 columnActivity_[i] = valueScaled*scaleR; 717 717 } 718 718 for (i=0;i<numberRows_;i++) { 719 double valueScaled = rowActivity_[i];719 CoinWorkDouble valueScaled = rowActivity_[i]; 720 720 rowActivity_[i] = valueScaled*scaleR; 721 721 } … … 762 762 } 763 763 int numberBad ; 764 double largestBound, smallestBound, minimumGap;765 double smallestObj, largestObj;764 CoinWorkDouble largestBound, smallestBound, minimumGap; 765 CoinWorkDouble smallestObj, largestObj; 766 766 int firstBad; 767 767 int modifiedBounds=0; … … 775 775 largestObj=0.0; 776 776 // If bounds are too close - fix 777 double fixTolerance = 1.1*primalTolerance();777 CoinWorkDouble fixTolerance = 1.1*primalTolerance(); 778 778 for (i=numberColumns_;i<numberColumns_+numberRows_;i++) { 779 double value;780 value = fabs(cost_[i]);779 CoinWorkDouble value; 780 value = CoinAbs(cost_[i]); 781 781 if (value>1.0e50) { 782 782 numberBad++; … … 805 805 } 806 806 if (lower_[i]>-1.0e100&&lower_[i]) { 807 value = fabs(lower_[i]);807 value = CoinAbs(lower_[i]); 808 808 if (value>largestBound) 809 809 largestBound=value; … … 812 812 } 813 813 if (upper_[i]<1.0e100&&upper_[i]) { 814 value = fabs(upper_[i]);814 value = CoinAbs(upper_[i]); 815 815 if (value>largestBound) 816 816 largestBound=value; … … 821 821 if (largestBound) 822 822 handler_->message(CLP_RIMSTATISTICS3,messages_) 823 <<s mallestBound824 << largestBound825 << minimumGap823 <<static_cast<double>(smallestBound) 824 <<static_cast<double>(largestBound) 825 <<static_cast<double>(minimumGap) 826 826 <<CoinMessageEol; 827 827 minimumGap=1.0e100; … … 829 829 largestBound=0.0; 830 830 for (i=0;i<numberColumns_;i++) { 831 double value;832 value = fabs(cost_[i]);831 CoinWorkDouble value; 832 value = CoinAbs(cost_[i]); 833 833 if (value>1.0e50) { 834 834 numberBad++; … … 857 857 } 858 858 if (lower_[i]>-1.0e100&&lower_[i]) { 859 value = fabs(lower_[i]);859 value = CoinAbs(lower_[i]); 860 860 if (value>largestBound) 861 861 largestBound=value; … … 864 864 } 865 865 if (upper_[i]<1.0e100&&upper_[i]) { 866 value = fabs(upper_[i]);866 value = CoinAbs(upper_[i]); 867 867 if (value>largestBound) 868 868 largestBound=value; … … 885 885 <<CoinMessageEol; 886 886 handler_->message(CLP_RIMSTATISTICS1,messages_) 887 <<s mallestObj888 << largestObj887 <<static_cast<double>(smallestObj) 888 <<static_cast<double>(largestObj) 889 889 <<CoinMessageEol; if (largestBound) 890 890 handler_->message(CLP_RIMSTATISTICS2,messages_) 891 <<s mallestBound892 << largestBound893 << minimumGap891 <<static_cast<double>(smallestBound) 892 <<static_cast<double>(largestBound) 893 <<static_cast<double>(minimumGap) 894 894 <<CoinMessageEol; 895 895 return true; … … 991 991 { 992 992 int iRow,iColumn; 993 CoinMemcpyN(cost_,numberColumns_,reducedCost_); 994 matrix_->transposeTimes(-1.0,dual_,reducedCost_); 993 CoinWorkDouble * reducedCost = reinterpret_cast<CoinWorkDouble *>(reducedCost_); 994 CoinWorkDouble * dual = reinterpret_cast<CoinWorkDouble *>(dual_); 995 CoinMemcpyN(cost_,numberColumns_,reducedCost); 996 matrix_->transposeTimes(-1.0,dual,reducedCost); 995 997 // Now modify reduced costs for quadratic 996 double quadraticOffset=quadraticDjs(reducedCost_,998 CoinWorkDouble quadraticOffset=quadraticDjs(reducedCost, 997 999 solution_,scaleFactor_); 998 1000 … … 1001 1003 sumPrimalInfeasibilities_=0.0; 1002 1004 sumDualInfeasibilities_=0.0; 1003 double dualTolerance = 10.0*dblParam_[ClpDualTolerance];1004 double primalTolerance = dblParam_[ClpPrimalTolerance];1005 double primalTolerance2 = 10.0*dblParam_[ClpPrimalTolerance];1005 CoinWorkDouble dualTolerance = 10.0*dblParam_[ClpDualTolerance]; 1006 CoinWorkDouble primalTolerance = dblParam_[ClpPrimalTolerance]; 1007 CoinWorkDouble primalTolerance2 = 10.0*dblParam_[ClpPrimalTolerance]; 1006 1008 worstComplementarity_=0.0; 1007 1009 complementarityGap_=0.0; … … 1009 1011 // Done scaled - use permanent regions for output 1010 1012 // but internal for bounds 1011 const double * lower = lower_+numberColumns_;1012 const double * upper = upper_+numberColumns_;1013 const CoinWorkDouble * lower = lower_+numberColumns_; 1014 const CoinWorkDouble * upper = upper_+numberColumns_; 1013 1015 for (iRow=0;iRow<numberRows_;iRow++) { 1014 double infeasibility=0.0;1015 double distanceUp = CoinMin(upper[iRow]-1016 rowActivity_[iRow],1.0e10);1017 double distanceDown = CoinMin(rowActivity_[iRow] -1018 lower[iRow],1.0e10);1016 CoinWorkDouble infeasibility=0.0; 1017 CoinWorkDouble distanceUp = CoinMin(upper[iRow]- 1018 rowActivity_[iRow],static_cast<CoinWorkDouble>(1.0e10)); 1019 CoinWorkDouble distanceDown = CoinMin(rowActivity_[iRow] - 1020 lower[iRow],static_cast<CoinWorkDouble>(1.0e10)); 1019 1021 if (distanceUp>primalTolerance2) { 1020 double value = dual_[iRow];1022 CoinWorkDouble value = dual[iRow]; 1021 1023 // should not be negative 1022 1024 if (value<-dualTolerance) { … … 1029 1031 } 1030 1032 if (distanceDown>primalTolerance2) { 1031 double value = dual_[iRow];1033 CoinWorkDouble value = dual[iRow]; 1032 1034 // should not be positive 1033 1035 if (value>dualTolerance) { … … 1051 1053 upper = upper_; 1052 1054 for (iColumn=0;iColumn<numberColumns_;iColumn++) { 1053 double infeasibility=0.0;1055 CoinWorkDouble infeasibility=0.0; 1054 1056 objectiveValue_ += cost_[iColumn]*columnActivity_[iColumn]; 1055 double distanceUp = CoinMin(upper[iColumn]-1056 columnActivity_[iColumn],1.0e10);1057 double distanceDown = CoinMin(columnActivity_[iColumn] -1058 lower[iColumn],1.0e10);1057 CoinWorkDouble distanceUp = CoinMin(upper[iColumn]- 1058 columnActivity_[iColumn],static_cast<CoinWorkDouble>(1.0e10)); 1059 CoinWorkDouble distanceDown = CoinMin(columnActivity_[iColumn] - 1060 lower[iColumn],static_cast<CoinWorkDouble>(1.0e10)); 1059 1061 if (distanceUp>primalTolerance2) { 1060 double value = reducedCost_[iColumn];1062 CoinWorkDouble value = reducedCost[iColumn]; 1061 1063 // should not be negative 1062 1064 if (value<-dualTolerance) { … … 1069 1071 } 1070 1072 if (distanceDown>primalTolerance2) { 1071 double value = reducedCost_[iColumn];1073 CoinWorkDouble value = reducedCost[iColumn]; 1072 1074 // should not be positive 1073 1075 if (value>dualTolerance) { … … 1088 1090 } 1089 1091 } 1092 #if COIN_LONG_WORK 1093 // ok as packs down 1094 CoinMemcpyN(reducedCost,numberColumns_,reducedCost_); 1095 #endif 1090 1096 // add in offset 1091 1097 objectiveValue_ += 0.5*quadraticOffset; … … 1142 1148 { 1143 1149 // Arrays for change in columns and rhs 1144 double * columnChange = new double[numberColumns_];1145 double * rowChange = new double[numberRows_];1150 CoinWorkDouble * columnChange = new CoinWorkDouble[numberColumns_]; 1151 CoinWorkDouble * rowChange = new CoinWorkDouble[numberRows_]; 1146 1152 CoinZeroN(columnChange,numberColumns_); 1147 1153 CoinZeroN(rowChange,numberRows_); 1148 1154 matrix_->times(1.0,columnChange,rowChange); 1149 1155 int i; 1150 double tolerance = primalTolerance();1156 CoinWorkDouble tolerance = primalTolerance(); 1151 1157 for (i=0;i<numberColumns_;i++) { 1152 1158 if (columnUpper_[i]<1.0e20||columnLower_[i]>-1.0e20) { … … 1154 1160 if (fixedOrFree(i)) { 1155 1161 if (columnActivity_[i]-columnLower_[i]<columnUpper_[i]-columnActivity_[i]) { 1156 double change = columnLower_[i]-columnActivity_[i];1157 if ( fabs(change)<tolerance) {1162 CoinWorkDouble change = columnLower_[i]-columnActivity_[i]; 1163 if (CoinAbs(change)<tolerance) { 1158 1164 if (reallyFix) 1159 1165 columnUpper_[i]=columnLower_[i]; … … 1162 1168 } 1163 1169 } else { 1164 double change = columnUpper_[i]-columnActivity_[i];1165 if ( fabs(change)<tolerance) {1170 CoinWorkDouble change = columnUpper_[i]-columnActivity_[i]; 1171 if (CoinAbs(change)<tolerance) { 1166 1172 if (reallyFix) 1167 1173 columnLower_[i]=columnUpper_[i]; … … 1177 1183 matrix_->times(1.0,columnChange,rowChange); 1178 1184 // If makes mess of things then don't do 1179 double newSum=0.0;1185 CoinWorkDouble newSum=0.0; 1180 1186 for (i=0;i<numberRows_;i++) { 1181 double value=rowActivity_[i]+rowChange[i];1187 CoinWorkDouble value=rowActivity_[i]+rowChange[i]; 1182 1188 if (value>rowUpper_[i]+tolerance) 1183 1189 newSum += value-rowUpper_[i]-tolerance; … … 1198 1204 if (fixedOrFree(i+numberColumns_)) { 1199 1205 if (rowActivity_[i]-rowLower_[i]<rowUpper_[i]-rowActivity_[i]) { 1200 double change = rowLower_[i]-rowActivity_[i];1201 if ( fabs(change)<tolerance) {1206 CoinWorkDouble change = rowLower_[i]-rowActivity_[i]; 1207 if (CoinAbs(change)<tolerance) { 1202 1208 if (reallyFix) 1203 1209 rowUpper_[i]=rowLower_[i]; … … 1205 1211 } 1206 1212 } else { 1207 double change = rowLower_[i]-rowActivity_[i];1208 if ( fabs(change)<tolerance) {1213 CoinWorkDouble change = rowLower_[i]-rowActivity_[i]; 1214 if (CoinAbs(change)<tolerance) { 1209 1215 if (reallyFix) 1210 1216 rowLower_[i]=rowUpper_[i]; … … 1223 1229 /* Modifies djs to allow for quadratic. 1224 1230 returns quadratic offset */ 1225 double1226 ClpInterior::quadraticDjs( double * djRegion, const double * solution,1227 double scaleFactor)1228 { 1229 double quadraticOffset=0.0;1231 CoinWorkDouble 1232 ClpInterior::quadraticDjs(CoinWorkDouble * djRegion, const CoinWorkDouble * solution, 1233 CoinWorkDouble scaleFactor) 1234 { 1235 CoinWorkDouble quadraticOffset=0.0; 1230 1236 #ifndef NO_RTTI 1231 1237 ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(objective_)); … … 1243 1249 int numberColumns = quadratic->getNumCols(); 1244 1250 for (int iColumn=0;iColumn<numberColumns;iColumn++) { 1245 double value = 0.0;1251 CoinWorkDouble value = 0.0; 1246 1252 for (CoinBigIndex j=columnQuadraticStart[iColumn]; 1247 1253 j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) { 1248 1254 int jColumn = columnQuadratic[j]; 1249 double valueJ = solution[jColumn];1250 double elementValue = quadraticElement[j];1255 CoinWorkDouble valueJ = solution[jColumn]; 1256 CoinWorkDouble elementValue = quadraticElement[j]; 1251 1257 //value += valueI*valueJ*elementValue; 1252 1258 value += valueJ*elementValue; -
trunk/Clp/src/ClpInterior.hpp
r1264 r1367 174 174 {algorithm_=value; } 175 175 /// Sum of dual infeasibilities 176 inline double sumDualInfeasibilities() const176 inline CoinWorkDouble sumDualInfeasibilities() const 177 177 { return sumDualInfeasibilities_;} 178 178 /// Sum of primal infeasibilities 179 inline double sumPrimalInfeasibilities() const179 inline CoinWorkDouble sumPrimalInfeasibilities() const 180 180 { return sumPrimalInfeasibilities_;} 181 181 /// dualObjective. 182 inline double dualObjective() const182 inline CoinWorkDouble dualObjective() const 183 183 { return dualObjective_;} 184 184 /// primalObjective. 185 inline double primalObjective() const185 inline CoinWorkDouble primalObjective() const 186 186 { return primalObjective_;} 187 187 /// diagonalNorm 188 inline double diagonalNorm() const188 inline CoinWorkDouble diagonalNorm() const 189 189 { return diagonalNorm_;} 190 190 /// linearPerturbation 191 inline double linearPerturbation() const191 inline CoinWorkDouble linearPerturbation() const 192 192 { return linearPerturbation_;} 193 inline void setLinearPerturbation( double value)193 inline void setLinearPerturbation(CoinWorkDouble value) 194 194 { linearPerturbation_=value;} 195 /// projectionTolerance 196 inline CoinWorkDouble projectionTolerance() const 197 { return projectionTolerance_;} 198 inline void setProjectionTolerance(CoinWorkDouble value) 199 { projectionTolerance_=value;} 195 200 /// diagonalPerturbation 196 inline double diagonalPerturbation() const201 inline CoinWorkDouble diagonalPerturbation() const 197 202 { return diagonalPerturbation_;} 198 inline void setDiagonalPerturbation( double value)203 inline void setDiagonalPerturbation(CoinWorkDouble value) 199 204 { diagonalPerturbation_=value;} 200 205 /// gamma 201 inline double gamma() const206 inline CoinWorkDouble gamma() const 202 207 { return gamma_;} 203 inline void setGamma( double value)208 inline void setGamma(CoinWorkDouble value) 204 209 { gamma_=value;} 205 210 /// delta 206 inline double delta() const211 inline CoinWorkDouble delta() const 207 212 { return delta_;} 208 inline void setDelta( double value)213 inline void setDelta(CoinWorkDouble value) 209 214 { delta_=value;} 210 215 /// ComplementarityGap 211 inline double complementarityGap() const216 inline CoinWorkDouble complementarityGap() const 212 217 { return complementarityGap_;} 213 218 //@} … … 216 221 //@{ 217 222 /// Largest error on Ax-b 218 inline double largestPrimalError() const223 inline CoinWorkDouble largestPrimalError() const 219 224 { return largestPrimalError_;} 220 225 /// Largest error on basic duals 221 inline double largestDualError() const226 inline CoinWorkDouble largestDualError() const 222 227 { return largestDualError_;} 223 228 /// Maximum iterations … … 234 239 void fixFixed(bool reallyFix=true); 235 240 /// Primal erturbation vector 236 inline double * primalR() const241 inline CoinWorkDouble * primalR() const 237 242 { return primalR_;} 238 243 /// Dual erturbation vector 239 inline double * dualR() const244 inline CoinWorkDouble * dualR() const 240 245 { return dualR_;} 241 246 //@} … … 260 265 //@{ 261 266 /// Raw objective value (so always minimize) 262 inline double rawObjectiveValue() const267 inline CoinWorkDouble rawObjectiveValue() const 263 268 { return objectiveValue_;} 264 269 /// Returns 1 if sequence indicates column … … 272 277 /** Modifies djs to allow for quadratic. 273 278 returns quadratic offset */ 274 double quadraticDjs(double * djRegion, const double * solution,275 double scaleFactor);279 CoinWorkDouble quadraticDjs(CoinWorkDouble * djRegion, const CoinWorkDouble * solution, 280 CoinWorkDouble scaleFactor); 276 281 277 282 /// To say a variable is fixed … … 370 375 //@{ 371 376 /// Largest error on Ax-b 372 double largestPrimalError_;377 CoinWorkDouble largestPrimalError_; 373 378 /// Largest error on basic duals 374 double largestDualError_;379 CoinWorkDouble largestDualError_; 375 380 /// Sum of dual infeasibilities 376 double sumDualInfeasibilities_;381 CoinWorkDouble sumDualInfeasibilities_; 377 382 /// Sum of primal infeasibilities 378 double sumPrimalInfeasibilities_;383 CoinWorkDouble sumPrimalInfeasibilities_; 379 384 /// Worst complementarity 380 double worstComplementarity_;385 CoinWorkDouble worstComplementarity_; 381 386 /// 382 387 public: 383 double xsize_;384 double zsize_;388 CoinWorkDouble xsize_; 389 CoinWorkDouble zsize_; 385 390 protected: 386 391 /// Working copy of lower bounds (Owner of arrays below) 387 double * lower_;392 CoinWorkDouble * lower_; 388 393 /// Row lower bounds - working copy 389 double * rowLowerWork_;394 CoinWorkDouble * rowLowerWork_; 390 395 /// Column lower bounds - working copy 391 double * columnLowerWork_;396 CoinWorkDouble * columnLowerWork_; 392 397 /// Working copy of upper bounds (Owner of arrays below) 393 double * upper_;398 CoinWorkDouble * upper_; 394 399 /// Row upper bounds - working copy 395 double * rowUpperWork_;400 CoinWorkDouble * rowUpperWork_; 396 401 /// Column upper bounds - working copy 397 double * columnUpperWork_;402 CoinWorkDouble * columnUpperWork_; 398 403 /// Working copy of objective 399 double * cost_;404 CoinWorkDouble * cost_; 400 405 public: 401 406 /// Rhs 402 double * rhs_;403 double * x_;404 double * y_;405 double * dj_;407 CoinWorkDouble * rhs_; 408 CoinWorkDouble * x_; 409 CoinWorkDouble * y_; 410 CoinWorkDouble * dj_; 406 411 protected: 407 412 /// Pointer to Lsqr object … … 411 416 /// Below here is standard barrier stuff 412 417 /// mu. 413 double mu_;418 CoinWorkDouble mu_; 414 419 /// objectiveNorm. 415 double objectiveNorm_;420 CoinWorkDouble objectiveNorm_; 416 421 /// rhsNorm. 417 double rhsNorm_;422 CoinWorkDouble rhsNorm_; 418 423 /// solutionNorm. 419 double solutionNorm_;424 CoinWorkDouble solutionNorm_; 420 425 /// dualObjective. 421 double dualObjective_;426 CoinWorkDouble dualObjective_; 422 427 /// primalObjective. 423 double primalObjective_;428 CoinWorkDouble primalObjective_; 424 429 /// diagonalNorm. 425 double diagonalNorm_;430 CoinWorkDouble diagonalNorm_; 426 431 /// stepLength 427 double stepLength_;432 CoinWorkDouble stepLength_; 428 433 /// linearPerturbation 429 double linearPerturbation_;434 CoinWorkDouble linearPerturbation_; 430 435 /// diagonalPerturbation 431 double diagonalPerturbation_;436 CoinWorkDouble diagonalPerturbation_; 432 437 // gamma from Saunders and Tomlin regularized 433 double gamma_;438 CoinWorkDouble gamma_; 434 439 // delta from Saunders and Tomlin regularized 435 double delta_;440 CoinWorkDouble delta_; 436 441 /// targetGap 437 double targetGap_;442 CoinWorkDouble targetGap_; 438 443 /// projectionTolerance 439 double projectionTolerance_;444 CoinWorkDouble projectionTolerance_; 440 445 /// maximumRHSError. maximum Ax 441 double maximumRHSError_;446 CoinWorkDouble maximumRHSError_; 442 447 /// maximumBoundInfeasibility. 443 double maximumBoundInfeasibility_;448 CoinWorkDouble maximumBoundInfeasibility_; 444 449 /// maximumDualError. 445 double maximumDualError_;450 CoinWorkDouble maximumDualError_; 446 451 /// diagonalScaleFactor. 447 double diagonalScaleFactor_;452 CoinWorkDouble diagonalScaleFactor_; 448 453 /// scaleFactor. For scaling objective 449 double scaleFactor_;454 CoinWorkDouble scaleFactor_; 450 455 /// actualPrimalStep 451 double actualPrimalStep_;456 CoinWorkDouble actualPrimalStep_; 452 457 /// actualDualStep 453 double actualDualStep_;458 CoinWorkDouble actualDualStep_; 454 459 /// smallestInfeasibility 455 double smallestInfeasibility_;460 CoinWorkDouble smallestInfeasibility_; 456 461 /// historyInfeasibility. 457 462 #define LENGTH_HISTORY 5 458 double historyInfeasibility_[LENGTH_HISTORY];463 CoinWorkDouble historyInfeasibility_[LENGTH_HISTORY]; 459 464 /// complementarityGap. 460 double complementarityGap_;465 CoinWorkDouble complementarityGap_; 461 466 /// baseObjectiveNorm 462 double baseObjectiveNorm_;467 CoinWorkDouble baseObjectiveNorm_; 463 468 /// worstDirectionAccuracy 464 double worstDirectionAccuracy_;469 CoinWorkDouble worstDirectionAccuracy_; 465 470 /// maximumRHSChange 466 double maximumRHSChange_;471 CoinWorkDouble maximumRHSChange_; 467 472 /// errorRegion. i.e. Ax 468 double * errorRegion_;473 CoinWorkDouble * errorRegion_; 469 474 /// rhsFixRegion. 470 double * rhsFixRegion_;475 CoinWorkDouble * rhsFixRegion_; 471 476 /// upperSlack 472 double * upperSlack_;477 CoinWorkDouble * upperSlack_; 473 478 /// lowerSlack 474 double * lowerSlack_;479 CoinWorkDouble * lowerSlack_; 475 480 /// diagonal 476 double * diagonal_;481 CoinWorkDouble * diagonal_; 477 482 /// solution 478 double * solution_;483 CoinWorkDouble * solution_; 479 484 /// work array 480 double * workArray_;485 CoinWorkDouble * workArray_; 481 486 /// delta X 482 double * deltaX_;487 CoinWorkDouble * deltaX_; 483 488 /// delta Y 484 double * deltaY_;489 CoinWorkDouble * deltaY_; 485 490 /// deltaZ. 486 double * deltaZ_;491 CoinWorkDouble * deltaZ_; 487 492 /// deltaW. 488 double * deltaW_;493 CoinWorkDouble * deltaW_; 489 494 /// deltaS. 490 double * deltaSU_;491 double * deltaSL_;495 CoinWorkDouble * deltaSU_; 496 CoinWorkDouble * deltaSL_; 492 497 /// Primal regularization array 493 double * primalR_;498 CoinWorkDouble * primalR_; 494 499 /// Dual regularization array 495 double * dualR_;500 CoinWorkDouble * dualR_; 496 501 /// rhs B 497 double * rhsB_;502 CoinWorkDouble * rhsB_; 498 503 /// rhsU. 499 double * rhsU_;504 CoinWorkDouble * rhsU_; 500 505 /// rhsL. 501 double * rhsL_;506 CoinWorkDouble * rhsL_; 502 507 /// rhsZ. 503 double * rhsZ_;508 CoinWorkDouble * rhsZ_; 504 509 /// rhsW. 505 double * rhsW_;510 CoinWorkDouble * rhsW_; 506 511 /// rhs C 507 double * rhsC_;512 CoinWorkDouble * rhsC_; 508 513 /// zVec 509 double * zVec_;514 CoinWorkDouble * zVec_; 510 515 /// wVec 511 double * wVec_;516 CoinWorkDouble * wVec_; 512 517 /// cholesky. 513 518 ClpCholeskyBase * cholesky_; -
trunk/Clp/src/ClpMatrixBase.cpp
r1034 r1367 605 605 abort(); 606 606 } 607 #if COIN_LONG_WORK 608 // For long double versions (aborts if not supported) 609 void 610 ClpMatrixBase::times(CoinWorkDouble scalar, 611 const CoinWorkDouble * x, CoinWorkDouble * y) const 612 { 613 std::cerr<<"long times not supported - ClpMatrixBase"<<std::endl; 614 abort(); 615 } 616 void 617 ClpMatrixBase::transposeTimes(CoinWorkDouble scalar, 618 const CoinWorkDouble * x, CoinWorkDouble * y) const 619 { 620 std::cerr<<"long transposeTimes not supported - ClpMatrixBase"<<std::endl; 621 abort(); 622 } 623 #endif -
trunk/Clp/src/ClpMatrixBase.hpp
r1302 r1367 5 5 6 6 #include "CoinPragma.hpp" 7 #include "CoinFinite.hpp" 7 8 8 9 #include "CoinPackedMatrix.hpp" … … 277 278 const double * columnScale, 278 279 double * spare=NULL) const; 280 #if COIN_LONG_WORK 281 // For long double versions (aborts if not supported) 282 virtual void times(CoinWorkDouble scalar, 283 const CoinWorkDouble * x, CoinWorkDouble * y) const ; 284 virtual void transposeTimes(CoinWorkDouble scalar, 285 const CoinWorkDouble * x, CoinWorkDouble * y) const ; 286 #endif 279 287 /** Return <code>x * scalar *A + y</code> in <code>z</code>. 280 288 Can use y as temporary array (will be empty at end) -
trunk/Clp/src/ClpMessage.cpp
r1034 r1367 82 82 {CLP_BARRIER_INFO,45,3,"Detail - %s"}, 83 83 {CLP_BARRIER_END,46,1,"At end primal/dual infeasibilities %g/%g, complementarity gap %g, objective %g"}, 84 {CLP_BARRIER_ACCURACY,47,2,"Relative error in phase %d, refinement %d is%g"},84 {CLP_BARRIER_ACCURACY,47,2,"Relative error in phase %d, %d refinements reduced from %g to %g"}, 85 85 {CLP_BARRIER_SAFE,48,2,"Initial safe primal value %g, objective norm %g"}, 86 86 {CLP_BARRIER_NEGATIVE_GAPS,49,3,"%d negative gaps summing to %g"}, -
trunk/Clp/src/ClpPackedMatrix.cpp
r1344 r1367 5961 5961 output->setPackedMode(true); 5962 5962 } 5963 #if COIN_LONG_WORK 5964 // For long double versions 5965 void 5966 ClpPackedMatrix::times(CoinWorkDouble scalar, 5967 const CoinWorkDouble * x, CoinWorkDouble * y) const 5968 { 5969 int iRow,iColumn; 5970 // get matrix data pointers 5971 const int * row = matrix_->getIndices(); 5972 const CoinBigIndex * columnStart = matrix_->getVectorStarts(); 5973 const double * elementByColumn = matrix_->getElements(); 5974 //memset(y,0,matrix_->getNumRows()*sizeof(double)); 5975 assert (((flags_&2)!=0)==matrix_->hasGaps()); 5976 if (!(flags_&2)) { 5977 for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) { 5978 CoinBigIndex j; 5979 CoinWorkDouble value = x[iColumn]; 5980 if (value) { 5981 CoinBigIndex start = columnStart[iColumn]; 5982 CoinBigIndex end = columnStart[iColumn+1]; 5983 value *= scalar; 5984 for (j=start; j<end; j++) { 5985 iRow=row[j]; 5986 y[iRow] += value*elementByColumn[j]; 5987 } 5988 } 5989 } 5990 } else { 5991 const int * columnLength = matrix_->getVectorLengths(); 5992 for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) { 5993 CoinBigIndex j; 5994 CoinWorkDouble value = x[iColumn]; 5995 if (value) { 5996 CoinBigIndex start = columnStart[iColumn]; 5997 CoinBigIndex end = start + columnLength[iColumn]; 5998 value *= scalar; 5999 for (j=start; j<end; j++) { 6000 iRow=row[j]; 6001 y[iRow] += value*elementByColumn[j]; 6002 } 6003 } 6004 } 6005 } 6006 } 6007 void 6008 ClpPackedMatrix::transposeTimes(CoinWorkDouble scalar, 6009 const CoinWorkDouble * x, CoinWorkDouble * y) const 6010 { 6011 int iColumn; 6012 // get matrix data pointers 6013 const int * row = matrix_->getIndices(); 6014 const CoinBigIndex * columnStart = matrix_->getVectorStarts(); 6015 const double * elementByColumn = matrix_->getElements(); 6016 if (!(flags_&2)) { 6017 if (scalar==-1.0) { 6018 CoinBigIndex start=columnStart[0]; 6019 for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) { 6020 CoinBigIndex j; 6021 CoinBigIndex next=columnStart[iColumn+1]; 6022 CoinWorkDouble value=y[iColumn]; 6023 for (j=start;j<next;j++) { 6024 int jRow=row[j]; 6025 value -= x[jRow]*elementByColumn[j]; 6026 } 6027 start=next; 6028 y[iColumn] = value; 6029 } 6030 } else { 6031 CoinBigIndex start=columnStart[0]; 6032 for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) { 6033 CoinBigIndex j; 6034 CoinBigIndex next=columnStart[iColumn+1]; 6035 CoinWorkDouble value=0.0; 6036 for (j=start;j<next;j++) { 6037 int jRow=row[j]; 6038 value += x[jRow]*elementByColumn[j]; 6039 } 6040 start=next; 6041 y[iColumn] += value*scalar; 6042 } 6043 } 6044 } else { 6045 const int * columnLength = matrix_->getVectorLengths(); 6046 for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) { 6047 CoinBigIndex j; 6048 CoinWorkDouble value=0.0; 6049 CoinBigIndex start = columnStart[iColumn]; 6050 CoinBigIndex end = start + columnLength[iColumn]; 6051 for (j=start; j<end; j++) { 6052 int jRow=row[j]; 6053 value += x[jRow]*elementByColumn[j]; 6054 } 6055 y[iColumn] += value*scalar; 6056 } 6057 } 6058 } 6059 #endif 5963 6060 #ifdef CLP_ALL_ONE_FILE 5964 6061 #undef reference -
trunk/Clp/src/ClpPackedMatrix.hpp
r1344 r1367 264 264 /// Sets up an effective RHS 265 265 void useEffectiveRhs(ClpSimplex * model); 266 #if COIN_LONG_WORK 267 // For long double versions 268 virtual void times(CoinWorkDouble scalar, 269 const CoinWorkDouble * x, CoinWorkDouble * y) const ; 270 virtual void transposeTimes(CoinWorkDouble scalar, 271 const CoinWorkDouble * x, CoinWorkDouble * y) const ; 272 #endif 266 273 //@} 267 274 -
trunk/Clp/src/ClpPredictorCorrector.cpp
r1366 r1367 33 33 fwrite(&numberRows,sizeof(int),1,fp); 34 34 fwrite(&numberColumns,sizeof(int),1,fp); 35 double dsave[20];35 CoinWorkDouble dsave[20]; 36 36 memset(dsave,0,sizeof(dsave)); 37 fwrite(dsave,sizeof( double),20,fp);37 fwrite(dsave,sizeof(CoinWorkDouble),20,fp); 38 38 int msave[20]; 39 39 memset(msave,0,sizeof(msave)); 40 40 msave[0]=numberIterations_; 41 41 fwrite(msave,sizeof(int),20,fp); 42 fwrite(dual_,sizeof( double),numberRows,fp);43 fwrite(errorRegion_,sizeof( double),numberRows,fp);44 fwrite(rhsFixRegion_,sizeof( double),numberRows,fp);45 fwrite(solution_,sizeof( double),numberColumns,fp);46 fwrite(solution_+numberColumns,sizeof( double),numberRows,fp);47 fwrite(diagonal_,sizeof( double),numberColumns,fp);48 fwrite(diagonal_+numberColumns,sizeof( double),numberRows,fp);49 fwrite(wVec_,sizeof( double),numberColumns,fp);50 fwrite(wVec_+numberColumns,sizeof( double),numberRows,fp);51 fwrite(zVec_,sizeof( double),numberColumns,fp);52 fwrite(zVec_+numberColumns,sizeof( double),numberRows,fp);53 fwrite(upperSlack_,sizeof( double),numberColumns,fp);54 fwrite(upperSlack_+numberColumns,sizeof( double),numberRows,fp);55 fwrite(lowerSlack_,sizeof( double),numberColumns,fp);56 fwrite(lowerSlack_+numberColumns,sizeof( double),numberRows,fp);42 fwrite(dual_,sizeof(CoinWorkDouble),numberRows,fp); 43 fwrite(errorRegion_,sizeof(CoinWorkDouble),numberRows,fp); 44 fwrite(rhsFixRegion_,sizeof(CoinWorkDouble),numberRows,fp); 45 fwrite(solution_,sizeof(CoinWorkDouble),numberColumns,fp); 46 fwrite(solution_+numberColumns,sizeof(CoinWorkDouble),numberRows,fp); 47 fwrite(diagonal_,sizeof(CoinWorkDouble),numberColumns,fp); 48 fwrite(diagonal_+numberColumns,sizeof(CoinWorkDouble),numberRows,fp); 49 fwrite(wVec_,sizeof(CoinWorkDouble),numberColumns,fp); 50 fwrite(wVec_+numberColumns,sizeof(CoinWorkDouble),numberRows,fp); 51 fwrite(zVec_,sizeof(CoinWorkDouble),numberColumns,fp); 52 fwrite(zVec_+numberColumns,sizeof(CoinWorkDouble),numberRows,fp); 53 fwrite(upperSlack_,sizeof(CoinWorkDouble),numberColumns,fp); 54 fwrite(upperSlack_+numberColumns,sizeof(CoinWorkDouble),numberRows,fp); 55 fwrite(lowerSlack_,sizeof(CoinWorkDouble),numberColumns,fp); 56 fwrite(lowerSlack_+numberColumns,sizeof(CoinWorkDouble),numberRows,fp); 57 57 fclose(fp); 58 58 } else { … … 61 61 } 62 62 #endif 63 static double eScale=1.0e27; 64 static double eBaseCaution=1.0e-12; 65 static double eBase=1.0e-12; 66 static double eRatio=1.0e40; 67 static double eRatioCaution=1.0e25; 68 static double eDiagonal=1.0e25; 69 static double eDiagonalCaution=1.0e18; 70 static double eExtra=1.0e-12; 63 // Could change on CLP_LONG_CHOLESKY or COIN_LONG_WORK? 64 static CoinWorkDouble eScale=1.0e27; 65 static CoinWorkDouble eBaseCaution=1.0e-12; 66 static CoinWorkDouble eBase=1.0e-12; 67 static CoinWorkDouble eRatio=1.0e40; 68 static CoinWorkDouble eRatioCaution=1.0e25; 69 static CoinWorkDouble eDiagonal=1.0e25; 70 static CoinWorkDouble eDiagonalCaution=1.0e18; 71 static CoinWorkDouble eExtra=1.0e-12; 71 72 72 73 // main function … … 81 82 return 2; 82 83 } 84 #if COIN_LONG_WORK 85 // reallocate some regions 86 double * dualSave = dual_; 87 dual_ = reinterpret_cast<double *>(new CoinWorkDouble[numberRows_]); 88 double * reducedCostSave = reducedCost_; 89 reducedCost_ = reinterpret_cast<double *>(new CoinWorkDouble[numberColumns_]); 90 #endif 83 91 //diagonalPerturbation_=1.0e-25; 84 92 ClpMatrixBase * saveMatrix = NULL; … … 165 173 return -1; 166 174 } 175 CoinWorkDouble * dualArray = reinterpret_cast<CoinWorkDouble *>(dual_); 167 176 // Could try centering steps without any original step i.e. just center 168 177 //firstFactorization(false); 169 CoinZeroN(dual _,numberRows_);178 CoinZeroN(dualArray,numberRows_); 170 179 multiplyAdd(solution_+numberColumns_,numberRows_,-1.0,errorRegion_,0.0); 171 180 matrix_->times(1.0,solution_,errorRegion_); 172 181 maximumRHSError_=maximumAbsElement(errorRegion_,numberRows_); 173 182 maximumBoundInfeasibility_=maximumRHSError_; 174 // double maximumDualError_=COIN_DBL_MAX;183 //CoinWorkDouble maximumDualError_=COIN_DBL_MAX; 175 184 //initialize 176 185 actualDualStep_=0.0; … … 184 193 int numberFixedTotal=numberFixed; 185 194 //int numberRows_DroppedBefore=0; 186 // double extra=eExtra;187 // double maximumPerturbation=COIN_DBL_MAX;195 //CoinWorkDouble extra=eExtra; 196 //CoinWorkDouble maximumPerturbation=COIN_DBL_MAX; 188 197 //constants for infeas interior point 189 const double beta2 = 0.99995;190 const double tau = 0.00002;191 double lastComplementarityGap=COIN_DBL_MAX;192 double lastStep=1.0;198 const CoinWorkDouble beta2 = 0.99995; 199 const CoinWorkDouble tau = 0.00002; 200 CoinWorkDouble lastComplementarityGap=COIN_DBL_MAX; 201 CoinWorkDouble lastStep=1.0; 193 202 // use to see if to take affine 194 double checkGap = COIN_DBL_MAX;203 CoinWorkDouble checkGap = COIN_DBL_MAX; 195 204 int lastGoodIteration=0; 196 double bestObjectiveGap=COIN_DBL_MAX; 197 double bestObjective=COIN_DBL_MAX; 205 CoinWorkDouble bestObjectiveGap=COIN_DBL_MAX; 206 CoinWorkDouble bestObjective=COIN_DBL_MAX; 207 int bestKilled=-1; 198 208 int saveIteration=-1; 209 int saveIteration2=-1; 199 210 bool sloppyOptimal=false; 200 double * savePi=NULL; 201 double * savePrimal=NULL; 211 CoinWorkDouble * savePi=NULL; 212 CoinWorkDouble * savePrimal=NULL; 213 CoinWorkDouble * savePi2=NULL; 214 CoinWorkDouble * savePrimal2=NULL; 202 215 // Extra regions for centering 203 double * saveX = new double[numberTotal];204 double * saveY = new double[numberRows_];205 double * saveZ = new double[numberTotal];206 double * saveW = new double[numberTotal];207 double * saveSL = new double[numberTotal];208 double * saveSU = new double[numberTotal];216 CoinWorkDouble * saveX = new CoinWorkDouble[numberTotal]; 217 CoinWorkDouble * saveY = new CoinWorkDouble[numberRows_]; 218 CoinWorkDouble * saveZ = new CoinWorkDouble[numberTotal]; 219 CoinWorkDouble * saveW = new CoinWorkDouble[numberTotal]; 220 CoinWorkDouble * saveSL = new CoinWorkDouble[numberTotal]; 221 CoinWorkDouble * saveSU = new CoinWorkDouble[numberTotal]; 209 222 // Save smallest mu used in primal dual moves 210 double smallestPrimalDualMu=COIN_DBL_MAX;211 double objScale = optimizationDirection_/223 CoinWorkDouble smallestPrimalDualMu=COIN_DBL_MAX; 224 CoinWorkDouble objScale = optimizationDirection_/ 212 225 (rhsScale_*objectiveScale_); 213 226 while (problemStatus_<0) { … … 231 244 handler_->message(CLP_BARRIER_ITERATION,messages_) 232 245 <<numberIterations_ 233 << primalObjective_*objScale- dblParam_[ClpObjOffset]234 << dualObjective_*objScale- dblParam_[ClpObjOffset]235 << complementarityGap_246 <<static_cast<double>(primalObjective_*objScale- dblParam_[ClpObjOffset]) 247 << static_cast<double>(dualObjective_*objScale- dblParam_[ClpObjOffset]) 248 <<static_cast<double>(complementarityGap_) 236 249 <<numberFixedTotal 237 250 <<cholesky_->rank() … … 253 266 //saveIteration=-1; 254 267 lastStep = CoinMin(actualPrimalStep_,actualDualStep_); 255 double goodGapChange;256 //#define KEEP_GOING_IF_FIXED 0268 CoinWorkDouble goodGapChange; 269 //#define KEEP_GOING_IF_FIXED 5 257 270 #ifndef KEEP_GOING_IF_FIXED 258 271 #define KEEP_GOING_IF_FIXED 10000 … … 265 278 goodGapChange=0.99; // make more likely to carry on 266 279 } 267 double gapO;268 double lastGood=bestObjectiveGap;280 CoinWorkDouble gapO; 281 CoinWorkDouble lastGood=bestObjectiveGap; 269 282 if (gonePrimalFeasible_&&goneDualFeasible_) { 270 double largestObjective;271 if ( fabs(primalObjective_)>fabs(dualObjective_)) {272 largestObjective = fabs(primalObjective_);283 CoinWorkDouble largestObjective; 284 if (CoinAbs(primalObjective_)>CoinAbs(dualObjective_)) { 285 largestObjective = CoinAbs(primalObjective_); 273 286 } else { 274 largestObjective = fabs(dualObjective_);287 largestObjective = CoinAbs(dualObjective_); 275 288 } 276 289 if (largestObjective<1.0) { 277 290 largestObjective=1.0; 278 291 } 279 gapO= fabs(primalObjective_-dualObjective_)/largestObjective;292 gapO=CoinAbs(primalObjective_-dualObjective_)/largestObjective; 280 293 handler_->message(CLP_BARRIER_OBJECTIVE_GAP,messages_) 281 << gapO294 <<static_cast<double>(gapO) 282 295 <<CoinMessageEol; 283 296 //start saving best 284 if (gapO<bestObjectiveGap) 297 bool saveIt=false; 298 if (gapO<bestObjectiveGap) { 285 299 bestObjectiveGap=gapO; 300 #ifndef SAVE_ON_OBJ 301 saveIt=true; 302 #endif 303 } 286 304 if (primalObjective_<bestObjective) { 305 bestObjective=primalObjective_; 306 #ifdef SAVE_ON_OBJ 307 saveIt=true; 308 #endif 309 } 310 if (numberFixedTotal>bestKilled) { 311 bestKilled=numberFixedTotal; 312 #if KEEP_GOING_IF_FIXED<10 313 saveIt=true; 314 #endif 315 } 316 if (saveIt) { 317 #if KEEP_GOING_IF_FIXED<10 318 printf("saving\n"); 319 #endif 287 320 saveIteration=numberIterations_; 288 bestObjective=primalObjective_;289 321 if (!savePi) { 290 savePi=new double[numberRows_];291 savePrimal = new double [numberTotal];292 } 293 CoinMemcpyN(dual _,numberRows_,savePi);322 savePi=new CoinWorkDouble[numberRows_]; 323 savePrimal = new CoinWorkDouble [numberTotal]; 324 } 325 CoinMemcpyN(dualArray,numberRows_,savePi); 294 326 CoinMemcpyN(solution_,numberTotal,savePrimal); 295 327 } else if(gapO>2.0*bestObjectiveGap) { … … 300 332 //std::cout <<"could stop"<<std::endl; 301 333 //gapO=0.0; 302 if ( fabs(primalObjective_-dualObjective_)<dualTolerance()) {334 if (CoinAbs(primalObjective_-dualObjective_)<dualTolerance()) { 303 335 gapO=0.0; 304 336 } … … 308 340 handler_->message(CLP_BARRIER_GONE_INFEASIBLE,messages_) 309 341 <<CoinMessageEol; 342 CoinWorkDouble scaledRHSError=maximumRHSError_/(solutionNorm_+10.0); 343 // save alternate 344 if (numberFixedTotal>bestKilled&& 345 maximumBoundInfeasibility_<1.0e-6&& 346 scaledRHSError<1.0e-2) { 347 bestKilled=numberFixedTotal; 348 #if KEEP_GOING_IF_FIXED<10 349 printf("saving alternate\n"); 350 #endif 351 saveIteration2=numberIterations_; 352 if (!savePi2) { 353 savePi2=new CoinWorkDouble[numberRows_]; 354 savePrimal2 = new CoinWorkDouble [numberTotal]; 355 } 356 CoinMemcpyN(dualArray,numberRows_,savePi2); 357 CoinMemcpyN(solution_,numberTotal,savePrimal2); 358 } 310 359 if (sloppyOptimal) { 311 360 // vaguely optimal 312 double scaledRHSError=maximumRHSError_/solutionNorm_;313 361 if (maximumBoundInfeasibility_>1.0e-2|| 314 362 scaledRHSError>1.0e-2|| … … 322 370 } else { 323 371 // not close to optimal but check if getting bad 324 double scaledRHSError=maximumRHSError_/solutionNorm_;372 CoinWorkDouble scaledRHSError=maximumRHSError_/(solutionNorm_+10.0); 325 373 if ((maximumBoundInfeasibility_>1.0e-1|| 326 374 scaledRHSError>1.0e-1|| … … 354 402 sloppyOptimal=true; 355 403 handler_->message(CLP_BARRIER_CLOSE_TO_OPTIMAL,messages_) 356 <<numberIterations_<< complementarityGap_404 <<numberIterations_<<static_cast<double>(complementarityGap_) 357 405 <<CoinMessageEol; 358 406 } … … 363 411 if (complementarityGap_>=1.05*lastComplementarityGap) { 364 412 handler_->message(CLP_BARRIER_COMPLEMENTARITY,messages_) 365 << complementarityGap_<<"increasing"413 <<static_cast<double>(complementarityGap_)<<"increasing" 366 414 <<CoinMessageEol; 367 415 if (saveIteration>=0&&sloppyOptimal) { … … 380 428 complementarityGap_<1.0e-3) { 381 429 handler_->message(CLP_BARRIER_COMPLEMENTARITY,messages_) 382 << complementarityGap_<<"not decreasing"430 <<static_cast<double>(complementarityGap_)<<"not decreasing" 383 431 <<CoinMessageEol; 384 432 if (gapO>0.75*lastGood&&numberFixed<KEEP_GOING_IF_FIXED) { … … 388 436 complementarityGap_<1.0e-6) { 389 437 handler_->message(CLP_BARRIER_COMPLEMENTARITY,messages_) 390 << complementarityGap_<<"not decreasing"438 <<static_cast<double>(complementarityGap_)<<"not decreasing" 391 439 <<CoinMessageEol; 392 440 break; … … 421 469 } 422 470 if (gapO<1.0e-9) { 423 double value=gapO*complementarityGap_;471 CoinWorkDouble value=gapO*complementarityGap_; 424 472 value*=actualPrimalStep_; 425 473 value*=actualDualStep_; … … 433 481 } 434 482 } 435 double nextGap=COIN_DBL_MAX;483 CoinWorkDouble nextGap=COIN_DBL_MAX; 436 484 int nextNumber=0; 437 485 int nextNumberItems=0; … … 441 489 //prepare for cholesky. Set up scaled diagonal in deltaX 442 490 // ** for efficiency may be better if scale factor known before 443 double norm2=0.0;444 double maximumValue;491 CoinWorkDouble norm2=0.0; 492 CoinWorkDouble maximumValue; 445 493 getNorms(diagonal_,numberTotal,maximumValue,norm2); 446 diagonalNorm_ = sqrt(norm2/numberComplementarityPairs_);494 diagonalNorm_ = CoinSqrt(norm2/numberComplementarityPairs_); 447 495 diagonalScaleFactor_=1.0; 448 double maximumAllowable=eScale;496 CoinWorkDouble maximumAllowable=eScale; 449 497 //scale so largest is less than allowable ? could do better 450 double factor=0.5;498 CoinWorkDouble factor=0.5; 451 499 while (maximumValue>maximumAllowable) { 452 500 diagonalScaleFactor_*=factor; … … 455 503 if (diagonalScaleFactor_!=1.0) { 456 504 handler_->message(CLP_BARRIER_SCALING,messages_) 457 <<"diagonal"<< diagonalScaleFactor_505 <<"diagonal"<<static_cast<double>(diagonalScaleFactor_) 458 506 <<CoinMessageEol; 459 507 diagonalNorm_*=diagonalScaleFactor_; … … 496 544 } 497 545 int phase=0; // predictor, corrector , primal dual 498 double directionAccuracy=0.0;546 CoinWorkDouble directionAccuracy=0.0; 499 547 bool doCorrector=true; 500 548 bool goodMove=true; … … 517 565 debugMove(0,1.0e-2,1.0e-2); 518 566 } 519 double affineGap=nextGap;567 CoinWorkDouble affineGap=nextGap; 520 568 int bestPhase=0; 521 double bestNextGap=nextGap;569 CoinWorkDouble bestNextGap=nextGap; 522 570 // ? 523 571 bestNextGap=CoinMax(nextGap,0.8*complementarityGap_); … … 526 574 if (complementarityGap_>1.0e-4*numberComplementarityPairs_) { 527 575 //std::cout <<"predicted duality gap "<<nextGap<<std::endl; 528 double part1=nextGap/numberComplementarityPairs_;576 CoinWorkDouble part1=nextGap/numberComplementarityPairs_; 529 577 part1=nextGap/numberComplementarityItems_; 530 double part2=nextGap/complementarityGap_;578 CoinWorkDouble part2=nextGap/complementarityGap_; 531 579 mu_=part1*part2*part2; 532 580 #if 0 533 double papermu =complementarityGap_/numberComplementarityPairs_;534 double affmu = nextGap/nextNumber;535 double sigma = pow(affmu/papermu,3);581 CoinWorkDouble papermu =complementarityGap_/numberComplementarityPairs_; 582 CoinWorkDouble affmu = nextGap/nextNumber; 583 CoinWorkDouble sigma = pow(affmu/papermu,3); 536 584 printf("mu %g, papermu %g, affmu %g, sigma %g sigmamu %g\n", 537 585 mu_,papermu,affmu,sigma,sigma*papermu); 538 586 #endif 539 587 //printf("paper mu %g\n",(nextGap*nextGap*nextGap)/(complementarityGap_*complementarityGap_* 540 // ( double) numberComplementarityPairs_));588 // (CoinWorkDouble) numberComplementarityPairs_)); 541 589 } else { 542 double phi;590 CoinWorkDouble phi; 543 591 if (numberComplementarityPairs_<=5000) { 544 phi=pow(static_cast< double> (numberComplementarityPairs_),2.0);592 phi=pow(static_cast<CoinWorkDouble> (numberComplementarityPairs_),2.0); 545 593 } else { 546 phi=pow(static_cast< double> (numberComplementarityPairs_),1.5);594 phi=pow(static_cast<CoinWorkDouble> (numberComplementarityPairs_),1.5); 547 595 if (phi<500.0*500.0) { 548 596 phi=500.0*500.0; … … 552 600 } 553 601 //save information 554 double product=affineProduct();602 CoinWorkDouble product=affineProduct(); 555 603 #if 0 556 604 //can we do corrector step? 557 double xx= complementarityGap_*(beta2-tau) +product;605 CoinWorkDouble xx= complementarityGap_*(beta2-tau) +product; 558 606 if (xx>0.0) { 559 double saveMu = mu_;560 double mu2=numberComplementarityPairs_;607 CoinWorkDouble saveMu = mu_; 608 CoinWorkDouble mu2=numberComplementarityPairs_; 561 609 mu2=xx/mu2; 562 610 if (mu2>mu_) { … … 593 641 bestNextGap=COIN_DBL_MAX; 594 642 } 595 // double floatNumber = 2.0*numberComplementarityPairs_;643 //CoinWorkDouble floatNumber = 2.0*numberComplementarityPairs_; 596 644 //floatNumber = 1.0*numberComplementarityItems_; 597 645 //mu_=nextGap/floatNumber; … … 615 663 CoinMemcpyN(deltaSU_,numberTotal,saveSU); 616 664 #ifdef HALVE 617 double savePrimalStep = actualPrimalStep_;618 double saveDualStep = actualDualStep_;619 double saveMu = mu_;665 CoinWorkDouble savePrimalStep = actualPrimalStep_; 666 CoinWorkDouble saveDualStep = actualDualStep_; 667 CoinWorkDouble saveMu = mu_; 620 668 #endif 621 669 //set up for next step 622 670 setupForSolve(phase); 623 double directionAccuracy2=findDirectionVector(phase);671 CoinWorkDouble directionAccuracy2=findDirectionVector(phase); 624 672 if (directionAccuracy2>worstDirectionAccuracy_) { 625 673 worstDirectionAccuracy_=directionAccuracy2; 626 674 } 627 double testValue=1.0e2*directionAccuracy;675 CoinWorkDouble testValue=1.0e2*directionAccuracy; 628 676 if (1.0e2*projectionTolerance_>testValue) { 629 677 testValue=1.0e2*projectionTolerance_; … … 644 692 if (goodMove) { 645 693 phase=1; 646 double norm = findStepLength(phase);694 CoinWorkDouble norm = findStepLength(phase); 647 695 nextGap = complementarityGap(nextNumber,nextNumberItems,1); 648 696 debugMove(1,actualPrimalStep_,actualDualStep_); … … 668 716 //printf("halve %d\n",nHalve); 669 717 nHalve++; 670 const double lambda=0.5;718 const CoinWorkDouble lambda=0.5; 671 719 for (i=0;i<numberRows_;i++) 672 720 deltaY_[i] = lambda*deltaY_[i]+(1.0-lambda)*saveY[i]; … … 699 747 if (!goodMove) { 700 748 // Just primal dual step 701 double floatNumber;749 CoinWorkDouble floatNumber; 702 750 floatNumber = 2.0*numberComplementarityPairs_; 703 751 //floatNumber = numberComplementarityItems_; 704 double saveMu=mu_; // use one from predictor corrector752 CoinWorkDouble saveMu=mu_; // use one from predictor corrector 705 753 mu_=complementarityGap_/floatNumber; 706 754 // If going well try small mu 707 mu_ *= sqrt((1.0-lastStep)/(1.0+10.0*lastStep));708 double mu1=mu_;709 double phi;755 mu_ *= CoinSqrt((1.0-lastStep)/(1.0+10.0*lastStep)); 756 CoinWorkDouble mu1=mu_; 757 CoinWorkDouble phi; 710 758 if (numberComplementarityPairs_<=500) { 711 phi=pow(static_cast< double> (numberComplementarityPairs_),2.0);759 phi=pow(static_cast<CoinWorkDouble> (numberComplementarityPairs_),2.0); 712 760 } else { 713 phi=pow(static_cast< double> (numberComplementarityPairs_),1.5);761 phi=pow(static_cast<CoinWorkDouble> (numberComplementarityPairs_),1.5); 714 762 if (phi<500.0*500.0) { 715 763 phi=500.0*500.0; … … 719 767 //printf("pd mu %g, alternate %g, smallest %g\n", 720 768 // mu_,mu1,smallestPrimalDualMu); 721 mu_ = sqrt(mu_*mu1);769 mu_ = CoinSqrt(mu_*mu1); 722 770 mu_=mu1; 723 771 if ((numberIterations_&1)==0||numberIterations_<10) … … 737 785 setupForSolve(2); 738 786 findDirectionVector(2); 739 double norm=findStepLength(2);787 CoinWorkDouble norm=findStepLength(2); 740 788 // just for debug 741 789 bestNextGap = complementarityGap_*1.0005; … … 760 808 smallestPrimalDualMu=mu_; 761 809 if (!goodMove) 762 mu_=nextGap / (static_cast< double> (nextNumber)*1.1);810 mu_=nextGap / (static_cast<CoinWorkDouble> (nextNumber)*1.1); 763 811 //if (quadraticObj) 764 812 //goodMove=true; … … 766 814 // Do centering steps 767 815 int numberTries=0; 768 double nextCenterGap=0.0;816 CoinWorkDouble nextCenterGap=0.0; 769 817 int numberGoodTries=0; 770 double originalDualStep=actualDualStep_;771 double originalPrimalStep=actualPrimalStep_;818 CoinWorkDouble originalDualStep=actualDualStep_; 819 CoinWorkDouble originalPrimalStep=actualPrimalStep_; 772 820 if (actualDualStep_>0.9&&actualPrimalStep_>0.9) 773 821 goodMove=false; // don't bother … … 781 829 CoinMemcpyN(deltaZ_,numberTotal,saveZ); 782 830 CoinMemcpyN(deltaW_,numberTotal,saveW); 783 double savePrimalStep = actualPrimalStep_;784 double saveDualStep = actualDualStep_;785 double saveMu = mu_;831 CoinWorkDouble savePrimalStep = actualPrimalStep_; 832 CoinWorkDouble saveDualStep = actualDualStep_; 833 CoinWorkDouble saveMu = mu_; 786 834 setupForSolve(3); 787 835 findDirectionVector(3); … … 789 837 debugMove(3,actualPrimalStep_,actualDualStep_); 790 838 //debugMove(3,1.0e-7,1.0e-7); 791 double xGap = complementarityGap(nextNumber,nextNumberItems,3);839 CoinWorkDouble xGap = complementarityGap(nextNumber,nextNumberItems,3); 792 840 // If one small then that's the one that counts 793 double checkDual=saveDualStep;794 double checkPrimal=savePrimalStep;841 CoinWorkDouble checkDual=saveDualStep; 842 CoinWorkDouble checkPrimal=savePrimalStep; 795 843 if (checkDual>5.0*checkPrimal) { 796 844 checkDual=2.0*checkPrimal; … … 848 896 delete [] saveSU; 849 897 if (savePi) { 850 std::cout<<"Restoring from iteration "<<saveIteration<<std::endl; 851 CoinMemcpyN(savePi,numberRows_,dual_); 852 CoinMemcpyN(savePrimal,numberTotal,solution_); 898 if (numberIterations_-saveIteration>20&& 899 numberIterations_-saveIteration2<5) { 900 #if KEEP_GOING_IF_FIXED<10 901 std::cout<<"Restoring2 from iteration "<<saveIteration2<<std::endl; 902 #endif 903 CoinMemcpyN(savePi2,numberRows_,dualArray); 904 CoinMemcpyN(savePrimal2,numberTotal,solution_); 905 } else { 906 #if KEEP_GOING_IF_FIXED<10 907 std::cout<<"Restoring from iteration "<<saveIteration<<std::endl; 908 #endif 909 CoinMemcpyN(savePi,numberRows_,dualArray); 910 CoinMemcpyN(savePrimal,numberTotal,solution_); 911 } 853 912 delete [] savePi; 854 913 delete [] savePrimal; 855 914 } 915 delete [] savePi2; 916 delete [] savePrimal2; 856 917 //recompute slacks 857 918 // Split out solution … … 861 922 //unscale objective 862 923 multiplyAdd(NULL,numberTotal,0.0,cost_,scaleFactor_); 863 multiplyAdd(NULL,numberRows_,0,dual _,scaleFactor_);924 multiplyAdd(NULL,numberRows_,0,dualArray,scaleFactor_); 864 925 checkSolution(); 865 CoinMemcpyN(reducedCost_,numberColumns_,dj_);926 //CoinMemcpyN(reducedCost_,numberColumns_,dj_); 866 927 // If quadratic use last solution 867 928 // Restore quadratic objective if necessary … … 872 933 } 873 934 handler_->message(CLP_BARRIER_END,messages_) 874 <<s umPrimalInfeasibilities_875 <<s umDualInfeasibilities_876 << complementarityGap_877 << objectiveValue()935 <<static_cast<double>(sumPrimalInfeasibilities_) 936 <<static_cast<double>(sumDualInfeasibilities_) 937 <<static_cast<double>(complementarityGap_) 938 <<static_cast<double>(objectiveValue()) 878 939 <<CoinMessageEol; 879 940 //#ifdef SOME_DEBUG … … 886 947 //std::cout<<"Primal objective "<<objectiveValue()<<std::endl; 887 948 //std::cout<<"maximum complementarity "<<worstComplementarity_<<std::endl; 949 #if COIN_LONG_WORK 950 // put back dual 951 delete [] dual_; 952 delete [] reducedCost_; 953 dual_ = dualSave; 954 reducedCost_=reducedCostSave; 955 #endif 888 956 //delete all temporary regions 889 957 deleteWorkingData(); 958 #if KEEP_GOING_IF_FIXED<10 959 #ifndef NDEBUG 960 { 961 static int kk=0; 962 char name[20]; 963 sprintf(name,"save.sol.%d",kk); 964 kk++; 965 printf("saving to file %s\n",name); 966 FILE * fp = fopen(name,"wb"); 967 int numberWritten; 968 numberWritten = fwrite(&numberColumns_,sizeof(int),1,fp); 969 assert (numberWritten==1); 970 numberWritten = fwrite(columnActivity_,sizeof(double),numberColumns_,fp); 971 assert (numberWritten==numberColumns_); 972 fclose(fp); 973 } 974 #endif 975 #endif 890 976 if (saveMatrix) { 891 977 // restore normal copy … … 899 985 // 1 corrector 900 986 // 2 primal dual 901 double ClpPredictorCorrector::findStepLength( int phase)987 CoinWorkDouble ClpPredictorCorrector::findStepLength( int phase) 902 988 { 903 double directionNorm=0.0;904 double maximumPrimalStep=COIN_DBL_MAX;905 double maximumDualStep=COIN_DBL_MAX;989 CoinWorkDouble directionNorm=0.0; 990 CoinWorkDouble maximumPrimalStep=COIN_DBL_MAX; 991 CoinWorkDouble maximumDualStep=COIN_DBL_MAX; 906 992 int numberTotal = numberRows_+numberColumns_; 907 double tolerance = 1.0e-12; 993 #if CLP_LONG_CHOLESKY<2 994 CoinWorkDouble tolerance = 1.0e-12; 995 #else 996 CoinWorkDouble tolerance = 1.0e-14; 997 #endif 908 998 int chosenPrimalSequence=-1; 909 999 int chosenDualSequence=-1; … … 911 1001 bool lowDual=false; 912 1002 // If done many iterations then allow to hit boundary 913 double hitTolerance;1003 CoinWorkDouble hitTolerance; 914 1004 //printf("objective norm %g\n",objectiveNorm_); 915 1005 if (numberIterations_<80||!gonePrimalFeasible_) … … 922 1012 for (iColumn=0;iColumn<numberTotal;iColumn++) { 923 1013 if (!flagged(iColumn)) { 924 double directionElement=deltaX_[iColumn]; 925 if (directionNorm<fabs(directionElement)) { 926 directionNorm=fabs(directionElement); 927 } 928 if (0) 929 printf("%d %g %g %g %g %g %g %g %g %g %g %g\n", 930 iColumn,solution_[iColumn],deltaX_[iColumn], 931 lowerSlack_[iColumn],deltaSL_[iColumn], 932 upperSlack_[iColumn],deltaSU_[iColumn], 933 dj_[iColumn], 934 zVec_[iColumn],deltaZ_[iColumn], 935 wVec_[iColumn],deltaW_[iColumn]); 1014 CoinWorkDouble directionElement=deltaX_[iColumn]; 1015 if (directionNorm<CoinAbs(directionElement)) { 1016 directionNorm=CoinAbs(directionElement); 1017 } 936 1018 if (lowerBound(iColumn)) { 937 double delta = - deltaSL_[iColumn];938 double z1 = deltaZ_[iColumn];939 double newZ = zVec_[iColumn]+z1;1019 CoinWorkDouble delta = - deltaSL_[iColumn]; 1020 CoinWorkDouble z1 = deltaZ_[iColumn]; 1021 CoinWorkDouble newZ = zVec_[iColumn]+z1; 940 1022 if (zVec_[iColumn]>tolerance) { 941 1023 if (zVec_[iColumn]<-z1*maximumDualStep) { … … 946 1028 } 947 1029 if (lowerSlack_[iColumn]<maximumPrimalStep*delta) { 948 double newStep=lowerSlack_[iColumn]/delta;1030 CoinWorkDouble newStep=lowerSlack_[iColumn]/delta; 949 1031 if (newStep>0.2||newZ<hitTolerance||delta>1.0e3||delta<=1.0e-6||dj_[iColumn]<hitTolerance) { 950 1032 maximumPrimalStep = newStep; … … 957 1039 } 958 1040 if (upperBound(iColumn)) { 959 double delta = - deltaSU_[iColumn];;960 double w1 = deltaW_[iColumn];961 double newT = wVec_[iColumn]+w1;1041 CoinWorkDouble delta = - deltaSU_[iColumn];; 1042 CoinWorkDouble w1 = deltaW_[iColumn]; 1043 CoinWorkDouble newT = wVec_[iColumn]+w1; 962 1044 if (wVec_[iColumn]>tolerance) { 963 1045 if (wVec_[iColumn]<-w1*maximumDualStep) { … … 968 1050 } 969 1051 if (upperSlack_[iColumn]<maximumPrimalStep*delta) { 970 double newStep=upperSlack_[iColumn]/delta;1052 CoinWorkDouble newStep=upperSlack_[iColumn]/delta; 971 1053 if (newStep>0.2||newT<hitTolerance||delta>1.0e3||delta<=1.0e-6||dj_[iColumn]>-hitTolerance) { 972 1054 maximumPrimalStep = newStep; … … 1026 1108 if (quadraticObj) { 1027 1109 // Use smaller unless very small 1028 double smallerStep=CoinMin(actualDualStep_,actualPrimalStep_);1110 CoinWorkDouble smallerStep=CoinMin(actualDualStep_,actualPrimalStep_); 1029 1111 if (smallerStep>0.0001) { 1030 1112 actualDualStep_=smallerStep; … … 1042 1124 // at first - just check better - if not 1043 1125 // Complementarity gap will be a*change*change + b*change +c 1044 double a=0.0;1045 double b=0.0;1046 double c=0.0;1126 CoinWorkDouble a=0.0; 1127 CoinWorkDouble b=0.0; 1128 CoinWorkDouble c=0.0; 1047 1129 /* SQUARE of dual infeasibility will be: 1048 1130 square of dj - ...... 1049 1131 */ 1050 double aq=0.0;1051 double bq=0.0;1052 double cq=0.0;1053 double gamma2 = gamma_*gamma_; // gamma*gamma will be added to diagonal1054 double * linearDjChange = new double[numberTotal];1132 CoinWorkDouble aq=0.0; 1133 CoinWorkDouble bq=0.0; 1134 CoinWorkDouble cq=0.0; 1135 CoinWorkDouble gamma2 = gamma_*gamma_; // gamma*gamma will be added to diagonal 1136 CoinWorkDouble * linearDjChange = new CoinWorkDouble[numberTotal]; 1055 1137 CoinZeroN(linearDjChange,numberColumns_); 1056 1138 multiplyAdd(deltaY_,numberRows_,1.0,linearDjChange+numberColumns_,0.0); … … 1060 1142 const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts(); 1061 1143 const int * columnQuadraticLength = quadratic->getVectorLengths(); 1062 double * quadraticElement = quadratic->getMutableElements();1144 CoinWorkDouble * quadraticElement = quadratic->getMutableElements(); 1063 1145 for (iColumn=0;iColumn<numberTotal;iColumn++) { 1064 double oldPrimal = solution_[iColumn];1146 CoinWorkDouble oldPrimal = solution_[iColumn]; 1065 1147 if (!flagged(iColumn)) { 1066 1148 if (lowerBound(iColumn)) { 1067 double change =oldPrimal+deltaX_[iColumn]-lowerSlack_[iColumn]-lower_[iColumn];1149 CoinWorkDouble change =oldPrimal+deltaX_[iColumn]-lowerSlack_[iColumn]-lower_[iColumn]; 1068 1150 c += lowerSlack_[iColumn]*zVec_[iColumn]; 1069 1151 b += lowerSlack_[iColumn]*deltaZ_[iColumn]+zVec_[iColumn]*change; … … 1071 1153 } 1072 1154 if (upperBound(iColumn)) { 1073 double change =upper_[iColumn]-oldPrimal-deltaX_[iColumn]-upperSlack_[iColumn];1155 CoinWorkDouble change =upper_[iColumn]-oldPrimal-deltaX_[iColumn]-upperSlack_[iColumn]; 1074 1156 c += upperSlack_[iColumn]*wVec_[iColumn]; 1075 1157 b += upperSlack_[iColumn]*deltaW_[iColumn]+wVec_[iColumn]*change; … … 1077 1159 } 1078 1160 // new djs are dj_ + change*value 1079 double djChange = linearDjChange[iColumn];1161 CoinWorkDouble djChange = linearDjChange[iColumn]; 1080 1162 if (iColumn<numberColumns_) { 1081 1163 for (CoinBigIndex j=columnQuadraticStart[iColumn]; 1082 1164 j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) { 1083 1165 int jColumn = columnQuadratic[j]; 1084 double changeJ = deltaX_[jColumn];1085 double elementValue = quadraticElement[j];1166 CoinWorkDouble changeJ = deltaX_[jColumn]; 1167 CoinWorkDouble elementValue = quadraticElement[j]; 1086 1168 djChange += changeJ*elementValue; 1087 1169 } 1088 1170 } 1089 double gammaTerm = gamma2;1171 CoinWorkDouble gammaTerm = gamma2; 1090 1172 if (primalR_) { 1091 1173 gammaTerm += primalR_[iColumn]; … … 1093 1175 djChange += gammaTerm; 1094 1176 // and dual infeasibility 1095 double oldInf = dj_[iColumn]-zVec_[iColumn]+wVec_[iColumn]+1177 CoinWorkDouble oldInf = dj_[iColumn]-zVec_[iColumn]+wVec_[iColumn]+ 1096 1178 gammaTerm*solution_[iColumn]; 1097 double changeInf = djChange-deltaZ_[iColumn]+deltaW_[iColumn];1179 CoinWorkDouble changeInf = djChange-deltaZ_[iColumn]+deltaW_[iColumn]; 1098 1180 cq += oldInf*oldInf; 1099 1181 bq += 2.0*oldInf*changeInf; … … 1108 1190 } 1109 1191 // new djs are dj_ + change*value 1110 double djChange = linearDjChange[iColumn];1192 CoinWorkDouble djChange = linearDjChange[iColumn]; 1111 1193 if (iColumn<numberColumns_) { 1112 1194 for (CoinBigIndex j=columnQuadraticStart[iColumn]; 1113 1195 j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) { 1114 1196 int jColumn = columnQuadratic[j]; 1115 double changeJ = deltaX_[jColumn];1116 double elementValue = quadraticElement[j];1197 CoinWorkDouble changeJ = deltaX_[jColumn]; 1198 CoinWorkDouble elementValue = quadraticElement[j]; 1117 1199 djChange += changeJ*elementValue; 1118 1200 } 1119 1201 } 1120 double gammaTerm = gamma2;1202 CoinWorkDouble gammaTerm = gamma2; 1121 1203 if (primalR_) { 1122 1204 gammaTerm += primalR_[iColumn]; … … 1124 1206 djChange += gammaTerm; 1125 1207 // and dual infeasibility 1126 double oldInf = dj_[iColumn]-zVec_[iColumn]+wVec_[iColumn]+1208 CoinWorkDouble oldInf = dj_[iColumn]-zVec_[iColumn]+wVec_[iColumn]+ 1127 1209 gammaTerm*solution_[iColumn]; 1128 double changeInf = djChange-deltaZ_[iColumn]+deltaW_[iColumn];1210 CoinWorkDouble changeInf = djChange-deltaZ_[iColumn]+deltaW_[iColumn]; 1129 1211 cq += oldInf*oldInf; 1130 1212 bq += 2.0*oldInf*changeInf; … … 1136 1218 // maybe use inf and do line search 1137 1219 // To check see if matches at current step 1138 double step=actualPrimalStep_;1139 //Current gap + solutionNorm_ * sqrt (sum square inf)1140 double multiplier = solutionNorm_;1220 CoinWorkDouble step=actualPrimalStep_; 1221 //Current gap + solutionNorm_ * CoinSqrt (sum square inf) 1222 CoinWorkDouble multiplier = solutionNorm_; 1141 1223 multiplier *= 0.01; 1142 1224 multiplier=1.0; 1143 double currentInf = multiplier*sqrt(cq);1144 double nextInf = multiplier*sqrt(CoinMax(cq+step*bq+step*step*aq,0.0));1145 double allowedIncrease=1.4;1225 CoinWorkDouble currentInf = multiplier*CoinSqrt(cq); 1226 CoinWorkDouble nextInf = multiplier*CoinSqrt(CoinMax(cq+step*bq+step*step*aq,0.0)); 1227 CoinWorkDouble allowedIncrease=1.4; 1146 1228 #ifdef SOME_DEBUG 1147 1229 printf("lin %g %g %g -> %g\n",a,b,c, … … 1161 1243 //cq = CoinMax(cq,10.0); 1162 1244 // convert to (x+q)*(x+q) = w 1163 double q = bq/(1.0*aq);1164 double w = CoinMax(q*q + (cq/aq)*(allowedIncrease-1.0),0.0);1165 w = sqrt(w);1166 double stepX = w-q;1245 CoinWorkDouble q = bq/(1.0*aq); 1246 CoinWorkDouble w = CoinMax(q*q + (cq/aq)*(allowedIncrease-1.0),0.0); 1247 w = CoinSqrt(w); 1248 CoinWorkDouble stepX = w-q; 1167 1249 step=stepX; 1168 1250 nextInf = 1169 multiplier* sqrt(CoinMax(cq+step*bq+step*step*aq,0.0));1251 multiplier*CoinSqrt(CoinMax(cq+step*bq+step*step*aq,0.0)); 1170 1252 #ifdef SOME_DEBUG 1171 1253 printf ("with step of %g dualInf is %g\n", … … 1180 1262 // minimize complementarity 1181 1263 // Complementarity gap will be a*change*change + b*change +c 1182 double a=0.0;1183 double b=0.0;1184 double c=0.0;1264 CoinWorkDouble a=0.0; 1265 CoinWorkDouble b=0.0; 1266 CoinWorkDouble c=0.0; 1185 1267 for (iColumn=0;iColumn<numberTotal;iColumn++) { 1186 double oldPrimal = solution_[iColumn];1268 CoinWorkDouble oldPrimal = solution_[iColumn]; 1187 1269 if (!flagged(iColumn)) { 1188 1270 if (lowerBound(iColumn)) { 1189 double change =oldPrimal+deltaX_[iColumn]-lowerSlack_[iColumn]-lower_[iColumn];1271 CoinWorkDouble change =oldPrimal+deltaX_[iColumn]-lowerSlack_[iColumn]-lower_[iColumn]; 1190 1272 c += lowerSlack_[iColumn]*zVec_[iColumn]; 1191 1273 b += lowerSlack_[iColumn]*deltaZ_[iColumn]+zVec_[iColumn]*change; … … 1193 1275 } 1194 1276 if (upperBound(iColumn)) { 1195 double change =upper_[iColumn]-oldPrimal-deltaX_[iColumn]-upperSlack_[iColumn];1277 CoinWorkDouble change =upper_[iColumn]-oldPrimal-deltaX_[iColumn]-upperSlack_[iColumn]; 1196 1278 c += upperSlack_[iColumn]*wVec_[iColumn]; 1197 1279 b += upperSlack_[iColumn]*deltaW_[iColumn]+wVec_[iColumn]*change; … … 1211 1293 // maybe use inf and do line search 1212 1294 // To check see if matches at current step 1213 double step=CoinMin(actualPrimalStep_,actualDualStep_);1214 double next = c+b*step+a*step*step;1295 CoinWorkDouble step=CoinMin(actualPrimalStep_,actualDualStep_); 1296 CoinWorkDouble next = c+b*step+a*step*step; 1215 1297 #ifdef SOME_DEBUG 1216 1298 printf("lin %g %g %g -> %g\n",a,b,c, … … 1230 1312 a=0.0; 1231 1313 b=0.0; 1232 double ratio = actualDualStep_/actualPrimalStep_;1314 CoinWorkDouble ratio = actualDualStep_/actualPrimalStep_; 1233 1315 for (iColumn=0;iColumn<numberTotal;iColumn++) { 1234 double oldPrimal = solution_[iColumn];1316 CoinWorkDouble oldPrimal = solution_[iColumn]; 1235 1317 if (!flagged(iColumn)) { 1236 1318 if (lowerBound(iColumn)) { 1237 double change =oldPrimal+deltaX_[iColumn]-lowerSlack_[iColumn]-lower_[iColumn];1319 CoinWorkDouble change =oldPrimal+deltaX_[iColumn]-lowerSlack_[iColumn]-lower_[iColumn]; 1238 1320 b += lowerSlack_[iColumn]*deltaZ_[iColumn]*ratio+zVec_[iColumn]*change; 1239 1321 a += deltaZ_[iColumn]*change*ratio; 1240 1322 } 1241 1323 if (upperBound(iColumn)) { 1242 double change =upper_[iColumn]-oldPrimal-deltaX_[iColumn]-upperSlack_[iColumn];1324 CoinWorkDouble change =upper_[iColumn]-oldPrimal-deltaX_[iColumn]-upperSlack_[iColumn]; 1243 1325 b += upperSlack_[iColumn]*deltaW_[iColumn]*ratio+wVec_[iColumn]*change; 1244 1326 a += deltaW_[iColumn]*change*ratio; … … 1250 1332 // To check see if matches at current step 1251 1333 step=actualPrimalStep_; 1252 double next2 = c+b*step+a*step*step;1334 CoinWorkDouble next2 = c+b*step+a*step*step; 1253 1335 if (next2>next) { 1254 1336 actualPrimalStep_=CoinMin(actualPrimalStep_,actualDualStep_); … … 1275 1357 #ifdef FULL_DEBUG 1276 1358 if (phase==3){ 1277 double minBeta = 0.1*mu_;1278 double maxBeta = 10.0*mu_;1359 CoinWorkDouble minBeta = 0.1*mu_; 1360 CoinWorkDouble maxBeta = 10.0*mu_; 1279 1361 for (iColumn=0;iColumn<numberRows_+numberColumns_;iColumn++) { 1280 1362 if (!flagged(iColumn)) { 1281 1363 if (lowerBound(iColumn)) { 1282 double change = -rhsL_[iColumn] + deltaX_[iColumn];1283 double dualValue=zVec_[iColumn]+actualDualStep_*deltaZ_[iColumn];1284 double primalValue=lowerSlack_[iColumn]+actualPrimalStep_*change;1285 double gapProduct=dualValue*primalValue;1364 CoinWorkDouble change = -rhsL_[iColumn] + deltaX_[iColumn]; 1365 CoinWorkDouble dualValue=zVec_[iColumn]+actualDualStep_*deltaZ_[iColumn]; 1366 CoinWorkDouble primalValue=lowerSlack_[iColumn]+actualPrimalStep_*change; 1367 CoinWorkDouble gapProduct=dualValue*primalValue; 1286 1368 if (delta2Z_[iColumn]<minBeta||delta2Z_[iColumn]>maxBeta) 1287 1369 printf("3lower %d primal %g, dual %g, gap %g, old gap %g\n", … … 1289 1371 } 1290 1372 if (upperBound(iColumn)) { 1291 double change = rhsU_[iColumn]-deltaX_[iColumn];1292 double dualValue=wVec_[iColumn]+actualDualStep_*deltaW_[iColumn];1293 double primalValue=upperSlack_[iColumn]+actualPrimalStep_*change;1294 double gapProduct=dualValue*primalValue;1373 CoinWorkDouble change = rhsU_[iColumn]-deltaX_[iColumn]; 1374 CoinWorkDouble dualValue=wVec_[iColumn]+actualDualStep_*deltaW_[iColumn]; 1375 CoinWorkDouble primalValue=upperSlack_[iColumn]+actualPrimalStep_*change; 1376 CoinWorkDouble gapProduct=dualValue*primalValue; 1295 1377 if (delta2W_[iColumn]<minBeta||delta2W_[iColumn]>maxBeta) 1296 1378 printf("3upper %d primal %g, dual %g, gap %g, old gap %g\n", … … 1303 1385 #ifdef SOME_DEBUG_not 1304 1386 { 1305 double largestL=0.0;1306 double smallestL=COIN_DBL_MAX;1307 double largestU=0.0;1308 double smallestU=COIN_DBL_MAX;1309 double sumL=0.0;1310 double sumU=0.0;1387 CoinWorkDouble largestL=0.0; 1388 CoinWorkDouble smallestL=COIN_DBL_MAX; 1389 CoinWorkDouble largestU=0.0; 1390 CoinWorkDouble smallestU=COIN_DBL_MAX; 1391 CoinWorkDouble sumL=0.0; 1392 CoinWorkDouble sumU=0.0; 1311 1393 int nL=0; 1312 1394 int nU=0; … … 1314 1396 if (!flagged(iColumn)) { 1315 1397 if (lowerBound(iColumn)) { 1316 double change = -rhsL_[iColumn] + deltaX_[iColumn];1317 double dualValue=zVec_[iColumn]+actualDualStep_*deltaZ_[iColumn];1318 double primalValue=lowerSlack_[iColumn]+actualPrimalStep_*change;1319 double gapProduct=dualValue*primalValue;1398 CoinWorkDouble change = -rhsL_[iColumn] + deltaX_[iColumn]; 1399 CoinWorkDouble dualValue=zVec_[iColumn]+actualDualStep_*deltaZ_[iColumn]; 1400 CoinWorkDouble primalValue=lowerSlack_[iColumn]+actualPrimalStep_*change; 1401 CoinWorkDouble gapProduct=dualValue*primalValue; 1320 1402 largestL = CoinMax(largestL,gapProduct); 1321 1403 smallestL = CoinMin(smallestL,gapProduct); … … 1324 1406 } 1325 1407 if (upperBound(iColumn)) { 1326 double change = rhsU_[iColumn]-deltaX_[iColumn];1327 double dualValue=wVec_[iColumn]+actualDualStep_*deltaW_[iColumn];1328 double primalValue=upperSlack_[iColumn]+actualPrimalStep_*change;1329 double gapProduct=dualValue*primalValue;1408 CoinWorkDouble change = rhsU_[iColumn]-deltaX_[iColumn]; 1409 CoinWorkDouble dualValue=wVec_[iColumn]+actualDualStep_*deltaW_[iColumn]; 1410 CoinWorkDouble primalValue=upperSlack_[iColumn]+actualPrimalStep_*change; 1411 CoinWorkDouble gapProduct=dualValue*primalValue; 1330 1412 largestU = CoinMax(largestU,gapProduct); 1331 1413 smallestU = CoinMin(smallestU,gapProduct); … … 1335 1417 } 1336 1418 } 1337 double mu = (sumL+sumU)/(static_cast<double> (nL+nU));1419 CoinWorkDouble mu = (sumL+sumU)/(static_cast<CoinWorkDouble> (nL+nU)); 1338 1420 1339 double minBeta = 0.1*mu;1340 double maxBeta = 10.0*mu;1421 CoinWorkDouble minBeta = 0.1*mu; 1422 CoinWorkDouble maxBeta = 10.0*mu; 1341 1423 int nBL=0; 1342 1424 int nAL=0; … … 1346 1428 if (!flagged(iColumn)) { 1347 1429 if (lowerBound(iColumn)) { 1348 double change = -rhsL_[iColumn] + deltaX_[iColumn];1349 double dualValue=zVec_[iColumn]+actualDualStep_*deltaZ_[iColumn];1350 double primalValue=lowerSlack_[iColumn]+actualPrimalStep_*change;1351 double gapProduct=dualValue*primalValue;1430 CoinWorkDouble change = -rhsL_[iColumn] + deltaX_[iColumn]; 1431 CoinWorkDouble dualValue=zVec_[iColumn]+actualDualStep_*deltaZ_[iColumn]; 1432 CoinWorkDouble primalValue=lowerSlack_[iColumn]+actualPrimalStep_*change; 1433 CoinWorkDouble gapProduct=dualValue*primalValue; 1352 1434 if (gapProduct<minBeta) 1353 1435 nBL++; … … 1359 1441 } 1360 1442 if (upperBound(iColumn)) { 1361 double change = rhsU_[iColumn]-deltaX_[iColumn];1362 double dualValue=wVec_[iColumn]+actualDualStep_*deltaW_[iColumn];1363 double primalValue=upperSlack_[iColumn]+actualPrimalStep_*change;1364 double gapProduct=dualValue*primalValue;1443 CoinWorkDouble change = rhsU_[iColumn]-deltaX_[iColumn]; 1444 CoinWorkDouble dualValue=wVec_[iColumn]+actualDualStep_*deltaW_[iColumn]; 1445 CoinWorkDouble primalValue=upperSlack_[iColumn]+actualPrimalStep_*change; 1446 CoinWorkDouble gapProduct=dualValue*primalValue; 1365 1447 if (gapProduct<minBeta) 1366 1448 nBU++; … … 1384 1466 /* Does solve. region1 is for deltaX (columns+rows), region2 for deltaPi (rows) */ 1385 1467 void 1386 ClpPredictorCorrector::solveSystem( double * region1, double * region2,1387 const double * region1In, const double * region2In,1388 const double * saveRegion1, const double * saveRegion2,1468 ClpPredictorCorrector::solveSystem(CoinWorkDouble * region1, CoinWorkDouble * region2, 1469 const CoinWorkDouble * region1In, const CoinWorkDouble * region2In, 1470 const CoinWorkDouble * saveRegion1, const CoinWorkDouble * saveRegion2, 1389 1471 bool gentleRefine) 1390 1472 { … … 1406 1488 multiplyAdd(region1+numberColumns_,numberRows_,-1.0,region2,1.0); 1407 1489 matrix_->times(1.0,region1,region2); 1408 double maximumRHS = maximumAbsElement(region2,numberRows_);1409 double scale=1.0;1410 double unscale=1.0;1490 CoinWorkDouble maximumRHS = maximumAbsElement(region2,numberRows_); 1491 CoinWorkDouble scale=1.0; 1492 CoinWorkDouble unscale=1.0; 1411 1493 if (maximumRHS>1.0e-30) { 1412 1494 if (maximumRHS<=0.5) { 1413 double factor=2.0;1495 CoinWorkDouble factor=2.0; 1414 1496 while (maximumRHS<=0.5) { 1415 1497 maximumRHS*=factor; … … 1417 1499 } /* endwhile */ 1418 1500 } else if (maximumRHS>=2.0&&maximumRHS<=COIN_DBL_MAX) { 1419 double factor=0.5;1501 CoinWorkDouble factor=0.5; 1420 1502 while (maximumRHS>=2.0) { 1421 1503 maximumRHS*=factor; … … 1444 1526 if (saveRegion2) { 1445 1527 //refine? 1446 double scaleX=1.0;1528 CoinWorkDouble scaleX=1.0; 1447 1529 if (gentleRefine) 1448 1530 scaleX=0.8; … … 1453 1535 } 1454 1536 // findDirectionVector. 1455 double ClpPredictorCorrector::findDirectionVector(const int phase)1537 CoinWorkDouble ClpPredictorCorrector::findDirectionVector(const int phase) 1456 1538 { 1457 double projectionTolerance=projectionTolerance_;1539 CoinWorkDouble projectionTolerance=projectionTolerance_; 1458 1540 //temporary 1459 1541 //projectionTolerance=1.0e-15; 1460 double errorCheck=0.9*maximumRHSError_/solutionNorm_;1542 CoinWorkDouble errorCheck=0.9*maximumRHSError_/solutionNorm_; 1461 1543 if (errorCheck>primalTolerance()) { 1462 1544 if (errorCheck<projectionTolerance) { … … 1468 1550 } 1469 1551 } 1470 double * newError = new double [numberRows_];1552 CoinWorkDouble * newError = new CoinWorkDouble [numberRows_]; 1471 1553 int numberTotal = numberRows_+numberColumns_; 1472 1554 //if flagged then entries zero so can do 1473 1555 // For KKT separate out 1474 double * region1Save=NULL;//for refinement1556 CoinWorkDouble * region1Save=NULL;//for refinement 1475 1557 int iColumn; 1476 1558 if (cholesky_->type()<20) { … … 1493 1575 } 1494 1576 bool goodSolve=false; 1495 double * regionSave=NULL;//for refinement1577 CoinWorkDouble * regionSave=NULL;//for refinement 1496 1578 int numberTries=0; 1497 double relativeError=COIN_DBL_MAX; 1498 double tryError=1.0e31; 1579 CoinWorkDouble relativeError=COIN_DBL_MAX; 1580 CoinWorkDouble tryError=1.0e31; 1581 CoinWorkDouble saveMaximum=0.0; 1582 double firstError=0.0; 1583 double lastError=0.0; 1499 1584 while (!goodSolve&&numberTries<30) { 1500 double lastError=relativeError;1585 CoinWorkDouble lastError=relativeError; 1501 1586 goodSolve=true; 1502 double maximumRHS; 1503 double saveMaximum; 1587 CoinWorkDouble maximumRHS; 1504 1588 maximumRHS = CoinMax(maximumAbsElement(deltaY_,numberRows_),1.0e-12); 1505 saveMaximum = maximumRHS; 1589 if (!numberTries) 1590 saveMaximum = maximumRHS; 1506 1591 if (cholesky_->type()<20) { 1507 1592 // no kkt 1508 double scale=1.0;1509 double unscale=1.0;1593 CoinWorkDouble scale=1.0; 1594 CoinWorkDouble unscale=1.0; 1510 1595 if (maximumRHS>1.0e-30) { 1511 1596 if (maximumRHS<=0.5) { 1512 double factor=2.0;1597 CoinWorkDouble factor=2.0; 1513 1598 while (maximumRHS<=0.5) { 1514 1599 maximumRHS*=factor; … … 1516 1601 } /* endwhile */ 1517 1602 } else if (maximumRHS>=2.0&&maximumRHS<=COIN_DBL_MAX) { 1518 double factor=0.5;1603 CoinWorkDouble factor=0.5; 1519 1604 while (maximumRHS>=2.0) { 1520 1605 maximumRHS*=factor; … … 1544 1629 if (numberTries) { 1545 1630 //refine? 1546 double scaleX=1.0;1631 CoinWorkDouble scaleX=1.0; 1547 1632 if (lastError>1.0e-5) 1548 1633 scaleX=0.8; … … 1567 1652 1568 1653 //now add in old Ax - doing extra checking 1569 double maximumRHSError=0.0;1570 double maximumRHSChange=0.0;1654 CoinWorkDouble maximumRHSError=0.0; 1655 CoinWorkDouble maximumRHSChange=0.0; 1571 1656 int iRow; 1572 1657 char * dropped = cholesky_->rowsDropped(); 1573 1658 for (iRow=0;iRow<numberRows_;iRow++) { 1574 1659 if (!dropped[iRow]) { 1575 double newValue=newError[iRow];1576 double oldValue=errorRegion_[iRow];1660 CoinWorkDouble newValue=newError[iRow]; 1661 CoinWorkDouble oldValue=errorRegion_[iRow]; 1577 1662 //severity of errors depend on signs 1578 1663 //**later */ 1579 if ( fabs(newValue)>maximumRHSChange) {1580 maximumRHSChange= fabs(newValue);1664 if (CoinAbs(newValue)>maximumRHSChange) { 1665 maximumRHSChange=CoinAbs(newValue); 1581 1666 } 1582 double result=newValue+oldValue;1583 if ( fabs(result)>maximumRHSError) {1584 maximumRHSError= fabs(result);1667 CoinWorkDouble result=newValue+oldValue; 1668 if (CoinAbs(result)>maximumRHSError) { 1669 maximumRHSError=CoinAbs(result); 1585 1670 } 1586 1671 newError[iRow]=result; 1587 1672 } else { 1588 double newValue=newError[iRow];1589 double oldValue=errorRegion_[iRow];1590 if ( fabs(newValue)>maximumRHSChange) {1591 maximumRHSChange= fabs(newValue);1673 CoinWorkDouble newValue=newError[iRow]; 1674 CoinWorkDouble oldValue=errorRegion_[iRow]; 1675 if (CoinAbs(newValue)>maximumRHSChange) { 1676 maximumRHSChange=CoinAbs(newValue); 1592 1677 } 1593 double result=newValue+oldValue;1678 CoinWorkDouble result=newValue+oldValue; 1594 1679 newError[iRow]=result; 1595 1680 //newError[iRow]=0.0; … … 1602 1687 if (relativeError>tryError) 1603 1688 relativeError=tryError; 1689 if (numberTries==1) 1690 firstError=relativeError; 1604 1691 if (relativeError<lastError) { 1692 lastError=relativeError; 1605 1693 maximumRHSChange_= maximumRHSChange; 1606 if (relativeError>1.0e-91607 ||numberTries>1) {1608 handler_->message(CLP_BARRIER_ACCURACY,messages_)1609 <<phase<<numberTries<<relativeError1610 <<CoinMessageEol;1611 }1612 1694 if (relativeError>projectionTolerance&&numberTries<=3) { 1613 1695 //try and refine … … 1617 1699 if (!goodSolve) { 1618 1700 if (!regionSave) { 1619 regionSave = new double [numberRows_];1701 regionSave = new CoinWorkDouble [numberRows_]; 1620 1702 if (cholesky_->type()>=20) 1621 region1Save = new double [numberTotal];1703 region1Save = new CoinWorkDouble [numberTotal]; 1622 1704 } 1623 1705 CoinMemcpyN(deltaY_,numberRows_,regionSave); … … 1653 1735 } else { 1654 1736 // disaster 1655 CoinFillN(deltaX_,numberTotal, 1.0);1656 CoinFillN(deltaY_,numberRows_, 1.0);1737 CoinFillN(deltaX_,numberTotal,static_cast<CoinWorkDouble>(1.0)); 1738 CoinFillN(deltaY_,numberRows_,static_cast<CoinWorkDouble>(1.0)); 1657 1739 printf("bad cholesky\n"); 1658 1740 } 1659 1741 } 1660 1742 } /* endwhile */ 1743 if (firstError>1.0e-9||numberTries>1) { 1744 handler_->message(CLP_BARRIER_ACCURACY,messages_) 1745 <<phase<<numberTries<<static_cast<double>(firstError) 1746 <<static_cast<double>(lastError) 1747 <<CoinMessageEol; 1748 } 1661 1749 delete [] regionSave; 1662 1750 delete [] region1Save; 1663 1751 delete [] newError; 1664 1752 // now rest 1665 double extra=eExtra;1753 CoinWorkDouble extra=eExtra; 1666 1754 //multiplyAdd(deltaY_,numberRows_,1.0,deltaW_+numberColumns_,0.0); 1667 1755 //CoinZeroN(deltaW_,numberColumns_); … … 1672 1760 deltaSL_[iColumn]=0.0; 1673 1761 deltaZ_[iColumn]=0.0; 1674 double dd=deltaW_[iColumn];1762 CoinWorkDouble dd=deltaW_[iColumn]; 1675 1763 deltaW_[iColumn]=0.0; 1676 1764 if (!flagged(iColumn)) { 1677 double deltaX = deltaX_[iColumn];1765 CoinWorkDouble deltaX = deltaX_[iColumn]; 1678 1766 if (lowerBound(iColumn)) { 1679 double zValue = rhsZ_[iColumn];1680 double gHat = zValue + zVec_[iColumn]*rhsL_[iColumn];1681 double slack = lowerSlack_[iColumn]+extra;1767 CoinWorkDouble zValue = rhsZ_[iColumn]; 1768 CoinWorkDouble gHat = zValue + zVec_[iColumn]*rhsL_[iColumn]; 1769 CoinWorkDouble slack = lowerSlack_[iColumn]+extra; 1682 1770 deltaSL_[iColumn] = -rhsL_[iColumn]+deltaX; 1683 1771 deltaZ_[iColumn]=(gHat-zVec_[iColumn]*deltaX)/slack; 1684 1772 } 1685 1773 if (upperBound(iColumn)) { 1686 double wValue = rhsW_[iColumn];1687 double hHat = wValue - wVec_[iColumn]*rhsU_[iColumn];1688 double slack = upperSlack_[iColumn]+extra;1774 CoinWorkDouble wValue = rhsW_[iColumn]; 1775 CoinWorkDouble hHat = wValue - wVec_[iColumn]*rhsU_[iColumn]; 1776 CoinWorkDouble slack = upperSlack_[iColumn]+extra; 1689 1777 deltaSU_[iColumn] = rhsU_[iColumn]-deltaX; 1690 1778 deltaW_[iColumn]=(hHat+wVec_[iColumn]*deltaX)/slack; … … 1692 1780 if (0) { 1693 1781 // different way of calculating 1694 double gamma2 = gamma_*gamma_;1695 double dZ=0.0;1696 double dW=0.0;1697 double zValue = rhsZ_[iColumn];1698 double gHat = zValue + zVec_[iColumn]*rhsL_[iColumn];1699 double slackL = lowerSlack_[iColumn]+extra;1700 double wValue = rhsW_[iColumn];1701 double hHat = wValue - wVec_[iColumn]*rhsU_[iColumn];1702 double slackU = upperSlack_[iColumn]+extra;1703 double q = rhsC_[iColumn]+gamma2 * deltaX +dd;1782 CoinWorkDouble gamma2 = gamma_*gamma_; 1783 CoinWorkDouble dZ=0.0; 1784 CoinWorkDouble dW=0.0; 1785 CoinWorkDouble zValue = rhsZ_[iColumn]; 1786 CoinWorkDouble gHat = zValue + zVec_[iColumn]*rhsL_[iColumn]; 1787 CoinWorkDouble slackL = lowerSlack_[iColumn]+extra; 1788 CoinWorkDouble wValue = rhsW_[iColumn]; 1789 CoinWorkDouble hHat = wValue - wVec_[iColumn]*rhsU_[iColumn]; 1790 CoinWorkDouble slackU = upperSlack_[iColumn]+extra; 1791 CoinWorkDouble q = rhsC_[iColumn]+gamma2 * deltaX +dd; 1704 1792 if (primalR_) 1705 1793 q += deltaX*primalR_[iColumn]; … … 1728 1816 } 1729 1817 #if 0 1730 double * check = new double[numberTotal];1818 CoinWorkDouble * check = new CoinWorkDouble[numberTotal]; 1731 1819 // Check out rhsC_ 1732 1820 multiplyAdd(deltaY_,numberRows_,-1.0,check+numberColumns_,0.0); … … 1736 1824 for (iColumn=0;iColumn<numberTotal;iColumn++) { 1737 1825 check[iColumn] += deltaZ_[iColumn]-deltaW_[iColumn]; 1738 if ( fabs(check[iColumn]-rhsC_[iColumn])>1.0e-3)1826 if (CoinAbs(check[iColumn]-rhsC_[iColumn])>1.0e-3) 1739 1827 printf("rhsC %d %g %g\n",iColumn,check[iColumn],rhsC_[iColumn]); 1740 1828 } … … 1743 1831 check[iColumn] += lowerSlack_[iColumn]*deltaZ_[iColumn]+ 1744 1832 zVec_[iColumn]*deltaSL_[iColumn]; 1745 if ( fabs(check[iColumn]-rhsZ_[iColumn])>1.0e-3)1833 if (CoinAbs(check[iColumn]-rhsZ_[iColumn])>1.0e-3) 1746 1834 printf("rhsZ %d %g %g\n",iColumn,check[iColumn],rhsZ_[iColumn]); 1747 1835 } … … 1750 1838 check[iColumn] += upperSlack_[iColumn]*deltaW_[iColumn]+ 1751 1839 wVec_[iColumn]*deltaSU_[iColumn]; 1752 if ( fabs(check[iColumn]-rhsW_[iColumn])>1.0e-3)1840 if (CoinAbs(check[iColumn]-rhsW_[iColumn])>1.0e-3) 1753 1841 printf("rhsW %d %g %g\n",iColumn,check[iColumn],rhsW_[iColumn]); 1754 1842 } … … 1762 1850 int numberTotal = numberRows_+numberColumns_; 1763 1851 int iColumn; 1764 double tolerance = primalTolerance();1852 CoinWorkDouble tolerance = primalTolerance(); 1765 1853 // See if quadratic objective 1766 1854 #ifndef NO_RTTI … … 1783 1871 clearFixed(iColumn); 1784 1872 } 1785 double initialValue =1.0e-12;1786 1873 1787 double maximumObjective=0.0;1788 double objectiveNorm2=0.0;1874 CoinWorkDouble maximumObjective=0.0; 1875 CoinWorkDouble objectiveNorm2=0.0; 1789 1876 getNorms(cost_,numberTotal,maximumObjective,objectiveNorm2); 1790 1877 if (!maximumObjective) { 1791 1878 maximumObjective=1.0; // objective all zero 1792 1879 } 1793 objectiveNorm2= sqrt(objectiveNorm2)/static_cast<double> (numberTotal);1880 objectiveNorm2=CoinSqrt(objectiveNorm2)/static_cast<CoinWorkDouble> (numberTotal); 1794 1881 objectiveNorm_=maximumObjective; 1795 1882 scaleFactor_=1.0; … … 1809 1896 if (quadraticObj) { 1810 1897 // If scaled then really scale matrix 1811 double scaleFactor =1898 CoinWorkDouble scaleFactor = 1812 1899 scaleFactor_*optimizationDirection_*objectiveScale_* 1813 1900 rhsScale_; … … 1819 1906 double * quadraticElement = quadratic->getMutableElements(); 1820 1907 int numberColumns = quadratic->getNumCols(); 1821 double scale = 1.0/scaleFactor;1908 CoinWorkDouble scale = 1.0/scaleFactor; 1822 1909 if (scalingFlag_>0&&rowScale_) { 1823 1910 for (int iColumn=0;iColumn<numberColumns;iColumn++) { 1824 double scaleI = columnScale_[iColumn]*scale;1911 CoinWorkDouble scaleI = columnScale_[iColumn]*scale; 1825 1912 for (CoinBigIndex j=columnQuadraticStart[iColumn]; 1826 1913 j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) { 1827 1914 int jColumn = columnQuadratic[j]; 1828 double scaleJ = columnScale_[jColumn];1915 CoinWorkDouble scaleJ = columnScale_[jColumn]; 1829 1916 quadraticElement[j] *= scaleI*scaleJ; 1830 objectiveNorm_ = CoinMax(objectiveNorm_, fabs(quadraticElement[j]));1917 objectiveNorm_ = CoinMax(objectiveNorm_,CoinAbs(quadraticElement[j])); 1831 1918 } 1832 1919 } … … 1837 1924 j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) { 1838 1925 quadraticElement[j] *= scale; 1839 objectiveNorm_ = CoinMax(objectiveNorm_, fabs(quadraticElement[j]));1926 objectiveNorm_ = CoinMax(objectiveNorm_,CoinAbs(quadraticElement[j])); 1840 1927 } 1841 1928 } … … 1848 1935 //DZ in lowerDual 1849 1936 //DW in upperDual 1850 double infiniteCheck=1.0e40;1851 // double fakeCheck=1.0e10;1937 CoinWorkDouble infiniteCheck=1.0e40; 1938 //CoinWorkDouble fakeCheck=1.0e10; 1852 1939 //use deltaX region for work region 1853 1940 for (iColumn=0;iColumn<numberTotal;iColumn++) { 1854 double primalValue = solution_[iColumn];1941 CoinWorkDouble primalValue = solution_[iColumn]; 1855 1942 clearFlagged(iColumn); 1856 1943 clearFixedOrFree(iColumn); … … 1863 1950 diagonal_[iColumn]=1.0; 1864 1951 deltaX_[iColumn]=1.0; 1865 double lowerValue=lower_[iColumn];1866 double upperValue=upper_[iColumn];1952 CoinWorkDouble lowerValue=lower_[iColumn]; 1953 CoinWorkDouble upperValue=upper_[iColumn]; 1867 1954 if (lowerValue>-infiniteCheck) { 1868 1955 if (upperValue<infiniteCheck) { … … 1953 2040 solveSystem(deltaX_,errorRegion_,solution_,NULL,NULL,NULL,false); 1954 2041 } 1955 initialValue=1.0e2;2042 CoinWorkDouble initialValue = 1.0e2; 1956 2043 if (rhsNorm_*1.0e-2>initialValue) { 1957 2044 initialValue=rhsNorm_*1.0e-2; 1958 2045 } 1959 2046 //initialValue = CoinMax(1.0,rhsNorm_); 1960 double smallestBoundDifference=COIN_DBL_MAX;1961 double * fakeSolution = deltaX_;2047 CoinWorkDouble smallestBoundDifference=COIN_DBL_MAX; 2048 CoinWorkDouble * fakeSolution = deltaX_; 1962 2049 for ( iColumn=0;iColumn<numberTotal;iColumn++) { 1963 2050 if (!flagged(iColumn)) { … … 1975 2062 solutionNorm_=1.0e-12; 1976 2063 handler_->message(CLP_BARRIER_SAFE,messages_) 1977 << initialValue<<objectiveNorm_2064 <<static_cast<double>(initialValue)<<static_cast<double>(objectiveNorm_) 1978 2065 <<CoinMessageEol; 1979 double extra=1.0e-10;1980 double largeGap=1.0e15;1981 // double safeObjectiveValue=2.0*objectiveNorm_;1982 double safeObjectiveValue=objectiveNorm_+1.0;1983 double safeFree=1.0e-1*initialValue;2066 CoinWorkDouble extra=1.0e-10; 2067 CoinWorkDouble largeGap=1.0e15; 2068 //CoinWorkDouble safeObjectiveValue=2.0*objectiveNorm_; 2069 CoinWorkDouble safeObjectiveValue=objectiveNorm_+1.0; 2070 CoinWorkDouble safeFree=1.0e-1*initialValue; 1984 2071 //printf("normal safe dual value of %g, primal value of %g\n", 1985 2072 // safeObjectiveValue,initialValue); … … 1988 2075 //printf("temp safe dual value of %g, primal value of %g\n", 1989 2076 // safeObjectiveValue,initialValue); 1990 double zwLarge=1.0e2*initialValue;2077 CoinWorkDouble zwLarge=1.0e2*initialValue; 1991 2078 //zwLarge=1.0e40; 1992 2079 if (cholesky_->choleskyCondition()<0.0&&cholesky_->type()<20) { … … 1996 2083 safeFree *= 10.0; 1997 2084 } 1998 double gamma2 = gamma_*gamma_; // gamma*gamma will be added to diagonal 1999 #if 0 2000 fakeSolution[0 ] = 0.072310129 ; 2001 fakeSolution[1 ] = 0.053083871; 2002 fakeSolution[2 ] = 0.178127; 2003 fakeSolution[3 ] = 0.13215151; 2004 fakeSolution[4 ] = 0.072715642; 2005 fakeSolution[5 ] = 0.15680727; 2006 fakeSolution[6 ] = 0.16841689; 2007 fakeSolution[7 ] = 0.093612798 ; 2008 fakeSolution[8 ] = 0.072774891 ; 2009 fakeSolution[9]=1.0; 2010 initialValue=1.0e-5; 2011 safeObjectiveValue=1.0e-5; 2012 #endif 2085 CoinWorkDouble gamma2 = gamma_*gamma_; // gamma*gamma will be added to diagonal 2013 2086 // First do primal side 2014 2087 for ( iColumn=0;iColumn<numberTotal;iColumn++) { 2015 2088 if (!flagged(iColumn)) { 2016 double lowerValue=lower_[iColumn];2017 double upperValue=upper_[iColumn];2018 double newValue;2019 double setPrimal=initialValue;2089 CoinWorkDouble lowerValue=lower_[iColumn]; 2090 CoinWorkDouble upperValue=upper_[iColumn]; 2091 CoinWorkDouble newValue; 2092 CoinWorkDouble setPrimal=initialValue; 2020 2093 if (quadraticObj) { 2021 2094 // perturb primal solution a bit … … 2026 2099 //upper and lower bounds 2027 2100 if (upperValue-lowerValue>2.0*setPrimal) { 2028 double fakeValue=fakeSolution[iColumn];2101 CoinWorkDouble fakeValue=fakeSolution[iColumn]; 2029 2102 if (fakeValue<lowerValue+setPrimal) { 2030 2103 fakeValue=lowerValue+setPrimal; … … 2039 2112 } else { 2040 2113 //just lower bound 2041 double fakeValue=fakeSolution[iColumn];2114 CoinWorkDouble fakeValue=fakeSolution[iColumn]; 2042 2115 if (fakeValue<lowerValue+setPrimal) { 2043 2116 fakeValue=lowerValue+setPrimal; … … 2048 2121 if (upperBound(iColumn)) { 2049 2122 //just upper bound 2050 double fakeValue=fakeSolution[iColumn];2123 CoinWorkDouble fakeValue=fakeSolution[iColumn]; 2051 2124 if (fakeValue>upperValue-setPrimal) { 2052 2125 fakeValue=upperValue-setPrimal; … … 2093 2166 quadraticElement = quadratic->getElements(); 2094 2167 } 2095 double quadraticNorm=0.0;2168 CoinWorkDouble quadraticNorm=0.0; 2096 2169 for ( iColumn=0;iColumn<numberTotal;iColumn++) { 2097 2170 if (!flagged(iColumn)) { 2098 double primalValue=solution_[iColumn];2099 double lowerValue=lower_[iColumn];2100 double upperValue=upper_[iColumn];2171 CoinWorkDouble primalValue=solution_[iColumn]; 2172 CoinWorkDouble lowerValue=lower_[iColumn]; 2173 CoinWorkDouble upperValue=upper_[iColumn]; 2101 2174 // Do dj 2102 double reducedCost=cost_[iColumn];2175 CoinWorkDouble reducedCost=cost_[iColumn]; 2103 2176 if (lowerBound(iColumn)) { 2104 2177 reducedCost+=linearPerturbation_; … … 2111 2184 j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) { 2112 2185 int jColumn = columnQuadratic[j]; 2113 double valueJ = solution_[jColumn];2114 double elementValue = quadraticElement[j];2186 CoinWorkDouble valueJ = solution_[jColumn]; 2187 CoinWorkDouble elementValue = quadraticElement[j]; 2115 2188 reducedCost += valueJ*elementValue; 2116 2189 } 2117 quadraticNorm = CoinMax(quadraticNorm, fabs(reducedCost));2190 quadraticNorm = CoinMax(quadraticNorm,CoinAbs(reducedCost)); 2118 2191 } 2119 2192 dj_[iColumn]=reducedCost; … … 2132 2205 for ( iColumn=0;iColumn<numberTotal;iColumn++) { 2133 2206 if (!flagged(iColumn)) { 2134 double primalValue=solution_[iColumn];2135 double lowerValue=lower_[iColumn];2136 double upperValue=upper_[iColumn];2137 double reducedCost=dj_[iColumn];2138 double low=0.0;2139 double high=0.0;2207 CoinWorkDouble primalValue=solution_[iColumn]; 2208 CoinWorkDouble lowerValue=lower_[iColumn]; 2209 CoinWorkDouble upperValue=upper_[iColumn]; 2210 CoinWorkDouble reducedCost=dj_[iColumn]; 2211 CoinWorkDouble low=0.0; 2212 CoinWorkDouble high=0.0; 2140 2213 if (lowerBound(iColumn)) { 2141 2214 if (upperBound(iColumn)) { … … 2148 2221 high=initialValue; 2149 2222 } 2150 double s = low+extra;2151 double ratioZ;2223 CoinWorkDouble s = low+extra; 2224 CoinWorkDouble ratioZ; 2152 2225 if (s<zwLarge) { 2153 2226 ratioZ=1.0; 2154 2227 } else { 2155 ratioZ= sqrt(zwLarge/s);2228 ratioZ=CoinSqrt(zwLarge/s); 2156 2229 } 2157 double t = high+extra;2158 double ratioT;2230 CoinWorkDouble t = high+extra; 2231 CoinWorkDouble ratioT; 2159 2232 if (t<zwLarge) { 2160 2233 ratioT=1.0; 2161 2234 } else { 2162 ratioT= sqrt(zwLarge/t);2235 ratioT=CoinSqrt(zwLarge/t); 2163 2236 } 2164 2237 //modify s and t … … 2179 2252 wVec_[iColumn]=CoinMax(-reducedCost , safeObjectiveValue*ratioT); 2180 2253 } 2181 double gammaTerm = gamma2;2254 CoinWorkDouble gammaTerm = gamma2; 2182 2255 if (primalR_) 2183 2256 gammaTerm += primalR_[iColumn]; … … 2188 2261 low=primalValue-lowerValue; 2189 2262 high=0.0; 2190 double s = low+extra;2191 double ratioZ;2263 CoinWorkDouble s = low+extra; 2264 CoinWorkDouble ratioZ; 2192 2265 if (s<zwLarge) { 2193 2266 ratioZ=1.0; 2194 2267 } else { 2195 ratioZ= sqrt(zwLarge/s);2268 ratioZ=CoinSqrt(zwLarge/s); 2196 2269 } 2197 2270 //modify s … … 2207 2280 wVec_[iColumn]=0.0; 2208 2281 } 2209 double gammaTerm = gamma2;2282 CoinWorkDouble gammaTerm = gamma2; 2210 2283 if (primalR_) 2211 2284 gammaTerm += primalR_[iColumn]; … … 2217 2290 low=0.0; 2218 2291 high=upperValue-primalValue; 2219 double t = high+extra;2220 double ratioT;2292 CoinWorkDouble t = high+extra; 2293 CoinWorkDouble ratioT; 2221 2294 if (t<zwLarge) { 2222 2295 ratioT=1.0; 2223 2296 } else { 2224 ratioT= sqrt(zwLarge/t);2297 ratioT=CoinSqrt(zwLarge/t); 2225 2298 } 2226 2299 //modify t … … 2236 2309 wVec_[iColumn]=CoinMax(-reducedCost , safeObjectiveValue*ratioT); 2237 2310 } 2238 double gammaTerm = gamma2;2311 CoinWorkDouble gammaTerm = gamma2; 2239 2312 if (primalR_) 2240 2313 gammaTerm += primalR_[iColumn]; … … 2249 2322 if (solution_[0]>0.0) { 2250 2323 for (int i=0;i<numberTotal;i++) 2251 printf("%d %.18g %.18g %.18g %.18g %.18g %.18g %.18g\n",i, fabs(solution_[i]),2252 diagonal_[i], fabs(dj_[i]),2324 printf("%d %.18g %.18g %.18g %.18g %.18g %.18g %.18g\n",i,CoinAbs(solution_[i]), 2325 diagonal_[i],CoinAbs(dj_[i]), 2253 2326 lowerSlack_[i],zVec_[i], 2254 2327 upperSlack_[i],wVec_[i]); 2255 2328 } else { 2256 2329 for (int i=0;i<numberTotal;i++) 2257 printf("%d %.18g %.18g %.18g %.18g %.18g %.18g %.18g\n",i, fabs(solution_[i]),2258 diagonal_[i], fabs(dj_[i]),2330 printf("%d %.18g %.18g %.18g %.18g %.18g %.18g %.18g\n",i,CoinAbs(solution_[i]), 2331 diagonal_[i],CoinAbs(dj_[i]), 2259 2332 upperSlack_[i],wVec_[i], 2260 2333 lowerSlack_[i],zVec_[i] ); … … 2266 2339 // complementarityGap. Computes gap 2267 2340 //phase 0=as is , 1 = after predictor , 2 after corrector 2268 double ClpPredictorCorrector::complementarityGap(int & numberComplementarityPairs,2341 CoinWorkDouble ClpPredictorCorrector::complementarityGap(int & numberComplementarityPairs, 2269 2342 int & numberComplementarityItems, 2270 2343 const int phase) 2271 2344 { 2272 double gap=0.0;2345 CoinWorkDouble gap=0.0; 2273 2346 //seems to be same coding for phase = 1 or 2 2274 2347 numberComplementarityPairs=0; 2275 2348 numberComplementarityItems=0; 2276 2349 int numberTotal = numberRows_+numberColumns_; 2277 double toleranceGap=0.0;2278 double largestGap=0.0;2279 double smallestGap=COIN_DBL_MAX;2350 CoinWorkDouble toleranceGap=0.0; 2351 CoinWorkDouble largestGap=0.0; 2352 CoinWorkDouble smallestGap=COIN_DBL_MAX; 2280 2353 //seems to be same coding for phase = 1 or 2 2281 2354 int numberNegativeGaps=0; 2282 double sumNegativeGap=0.0;2283 double largeGap=1.0e2*solutionNorm_;2355 CoinWorkDouble sumNegativeGap=0.0; 2356 CoinWorkDouble largeGap=1.0e2*solutionNorm_; 2284 2357 if (largeGap<1.0e10) { 2285 2358 largeGap=1.0e10; 2286 2359 } 2287 2360 largeGap=1.0e30; 2288 double dualTolerance = dblParam_[ClpDualTolerance];2289 double primalTolerance = dblParam_[ClpPrimalTolerance];2361 CoinWorkDouble dualTolerance = dblParam_[ClpDualTolerance]; 2362 CoinWorkDouble primalTolerance = dblParam_[ClpPrimalTolerance]; 2290 2363 dualTolerance=dualTolerance/scaleFactor_; 2291 2364 for (int iColumn=0;iColumn<numberTotal;iColumn++) { … … 2293 2366 numberComplementarityPairs++; 2294 2367 //can collapse as if no lower bound both zVec and deltaZ 0.0 2295 double newZ=0.0;2296 double newW=0.0;2368 CoinWorkDouble newZ=0.0; 2369 CoinWorkDouble newW=0.0; 2297 2370 if (lowerBound(iColumn)) { 2298 2371 numberComplementarityItems++; 2299 double dualValue;2300 double primalValue;2372 CoinWorkDouble dualValue; 2373 CoinWorkDouble primalValue; 2301 2374 if (!phase) { 2302 2375 dualValue=zVec_[iColumn]; 2303 2376 primalValue=lowerSlack_[iColumn]; 2304 2377 } else { 2305 double change;2378 CoinWorkDouble change; 2306 2379 change =solution_[iColumn]+deltaX_[iColumn]-lowerSlack_[iColumn]-lower_[iColumn]; 2307 2380 dualValue=zVec_[iColumn]+actualDualStep_*deltaZ_[iColumn]; … … 2313 2386 primalValue=largeGap; 2314 2387 } 2315 double gapProduct=dualValue*primalValue;2388 CoinWorkDouble gapProduct=dualValue*primalValue; 2316 2389 if (gapProduct<0.0) { 2317 2390 //cout<<"negative gap component "<<iColumn<<" "<<dualValue<<" "<< … … 2334 2407 if (upperBound(iColumn)) { 2335 2408 numberComplementarityItems++; 2336 double dualValue;2337 double primalValue;2409 CoinWorkDouble dualValue; 2410 CoinWorkDouble primalValue; 2338 2411 if (!phase) { 2339 2412 dualValue=wVec_[iColumn]; 2340 2413 primalValue=upperSlack_[iColumn]; 2341 2414 } else { 2342 double change;2415 CoinWorkDouble change; 2343 2416 change =upper_[iColumn]-solution_[iColumn]-deltaX_[iColumn]-upperSlack_[iColumn]; 2344 2417 dualValue=wVec_[iColumn]+actualDualStep_*deltaW_[iColumn]; … … 2350 2423 primalValue=largeGap; 2351 2424 } 2352 double gapProduct=dualValue*primalValue;2425 CoinWorkDouble gapProduct=dualValue*primalValue; 2353 2426 if (gapProduct<0.0) { 2354 2427 //cout<<"negative gap component "<<iColumn<<" "<<dualValue<<" "<< … … 2374 2447 if (!phase&&numberNegativeGaps) { 2375 2448 handler_->message(CLP_BARRIER_NEGATIVE_GAPS,messages_) 2376 <<numberNegativeGaps<<sumNegativeGap 2449 <<numberNegativeGaps<<static_cast<double>(sumNegativeGap) 2377 2450 <<CoinMessageEol; 2378 2451 } … … 2393 2466 void ClpPredictorCorrector::setupForSolve(const int phase) 2394 2467 { 2395 double extra =eExtra;2468 CoinWorkDouble extra =eExtra; 2396 2469 int numberTotal = numberRows_ + numberColumns_; 2397 2470 int iColumn; … … 2399 2472 printf("phase %d in setupForSolve, mu %.18g\n",phase,mu_); 2400 2473 #endif 2401 double gamma2 = gamma_*gamma_; // gamma*gamma will be added to diagonal 2474 CoinWorkDouble gamma2 = gamma_*gamma_; // gamma*gamma will be added to diagonal 2475 CoinWorkDouble * dualArray = reinterpret_cast<CoinWorkDouble *>(dual_); 2402 2476 switch (phase) { 2403 2477 case 0: … … 2405 2479 if (delta_||dualR_) { 2406 2480 // add in regularization 2407 double delta2 = delta_*delta_;2481 CoinWorkDouble delta2 = delta_*delta_; 2408 2482 for (int iRow=0;iRow<numberRows_;iRow++) { 2409 rhsB_[iRow] -= delta2*dual _[iRow];2483 rhsB_[iRow] -= delta2*dualArray[iRow]; 2410 2484 if (dualR_) 2411 rhsB_[iRow] -= dualR_[iRow]*dual _[iRow];2485 rhsB_[iRow] -= dualR_[iRow]*dualArray[iRow]; 2412 2486 } 2413 2487 } … … 2445 2519 #if 0 2446 2520 for (int i=0;i<3;i++) { 2447 if (! fabs(rhsZ_[i]))2521 if (!CoinAbs(rhsZ_[i])) 2448 2522 rhsZ_[i]=0.0; 2449 if (! fabs(rhsW_[i]))2523 if (!CoinAbs(rhsW_[i])) 2450 2524 rhsW_[i]=0.0; 2451 if (! fabs(rhsU_[i]))2525 if (!CoinAbs(rhsU_[i])) 2452 2526 rhsU_[i]=0.0; 2453 if (! fabs(rhsL_[i]))2527 if (!CoinAbs(rhsL_[i])) 2454 2528 rhsL_[i]=0.0; 2455 2529 } … … 2499 2573 #if 0 2500 2574 for (int i=0;i<numberTotal;i++) { 2501 if (! fabs(rhsZ_[i]))2575 if (!CoinAbs(rhsZ_[i])) 2502 2576 rhsZ_[i]=0.0; 2503 if (! fabs(rhsW_[i]))2577 if (!CoinAbs(rhsW_[i])) 2504 2578 rhsW_[i]=0.0; 2505 if (! fabs(rhsU_[i]))2579 if (!CoinAbs(rhsU_[i])) 2506 2580 rhsU_[i]=0.0; 2507 if (! fabs(rhsL_[i]))2581 if (!CoinAbs(rhsL_[i])) 2508 2582 rhsL_[i]=0.0; 2509 2583 } 2510 2584 if (solution_[0]>0.0) { 2511 2585 for (int i=0;i<numberTotal;i++) 2512 printf("%d %.18g %.18g %.18g %.18g %.18g %.18g %.18g\n",i, fabs(solution_[i]),2513 diagonal_[i], fabs(dj_[i]),2586 printf("%d %.18g %.18g %.18g %.18g %.18g %.18g %.18g\n",i,CoinAbs(solution_[i]), 2587 diagonal_[i],CoinAbs(dj_[i]), 2514 2588 lowerSlack_[i],zVec_[i], 2515 2589 upperSlack_[i],wVec_[i]); 2516 2590 for (int i=0;i<numberTotal;i++) 2517 printf("%d %.18g %.18g %.18g %.18g %.18g\n",i, fabs(rhsC_[i]),2591 printf("%d %.18g %.18g %.18g %.18g %.18g\n",i,CoinAbs(rhsC_[i]), 2518 2592 rhsZ_[i],rhsL_[i], 2519 2593 rhsW_[i],rhsU_[i]); 2520 2594 } else { 2521 2595 for (int i=0;i<numberTotal;i++) 2522 printf("%d %.18g %.18g %.18g %.18g %.18g %.18g %.18g\n",i, fabs(solution_[i]),2523 diagonal_[i], fabs(dj_[i]),2596 printf("%d %.18g %.18g %.18g %.18g %.18g %.18g %.18g\n",i,CoinAbs(solution_[i]), 2597 diagonal_[i],CoinAbs(dj_[i]), 2524 2598 upperSlack_[i],wVec_[i], 2525 2599 lowerSlack_[i],zVec_[i] ); 2526 2600 for (int i=0;i<numberTotal;i++) 2527 printf("%d %.18g %.18g %.18g %.18g %.18g\n",i, fabs(rhsC_[i]),2601 printf("%d %.18g %.18g %.18g %.18g %.18g\n",i,CoinAbs(rhsC_[i]), 2528 2602 rhsW_[i],rhsU_[i], 2529 2603 rhsZ_[i],rhsL_[i]); … … 2549 2623 case 3: 2550 2624 { 2551 double minBeta = 0.1*mu_;2552 double maxBeta = 10.0*mu_;2553 double dualStep = CoinMin(1.0,actualDualStep_+0.1);2554 double primalStep = CoinMin(1.0,actualPrimalStep_+0.1);2625 CoinWorkDouble minBeta = 0.1*mu_; 2626 CoinWorkDouble maxBeta = 10.0*mu_; 2627 CoinWorkDouble dualStep = CoinMin(1.0,actualDualStep_+0.1); 2628 CoinWorkDouble primalStep = CoinMin(1.0,actualPrimalStep_+0.1); 2555 2629 #ifdef SOME_DEBUG 2556 2630 printf("good complementarity range %g to %g\n",minBeta,maxBeta); … … 2561 2635 if (!flagged(iColumn)) { 2562 2636 if (lowerBound(iColumn)) { 2563 double change = -rhsL_[iColumn] + deltaX_[iColumn];2564 double dualValue=zVec_[iColumn]+dualStep*deltaZ_[iColumn];2565 double primalValue=lowerSlack_[iColumn]+primalStep*change;2566 double gapProduct=dualValue*primalValue;2637 CoinWorkDouble change = -rhsL_[iColumn] + deltaX_[iColumn]; 2638 CoinWorkDouble dualValue=zVec_[iColumn]+dualStep*deltaZ_[iColumn]; 2639 CoinWorkDouble primalValue=lowerSlack_[iColumn]+primalStep*change; 2640 CoinWorkDouble gapProduct=dualValue*primalValue; 2567 2641 if (gapProduct>0.0&&dualValue<0.0) 2568 2642 gapProduct = - gapProduct; … … 2573 2647 iColumn,primalValue,dualValue,gapProduct); 2574 2648 #endif 2575 double value= 0.0;2649 CoinWorkDouble value= 0.0; 2576 2650 if (gapProduct<minBeta) { 2577 2651 value= 2.0*(minBeta-gapProduct); … … 2586 2660 } 2587 2661 if (upperBound(iColumn)) { 2588 double change = rhsU_[iColumn]-deltaX_[iColumn];2589 double dualValue=wVec_[iColumn]+dualStep*deltaW_[iColumn];2590 double primalValue=upperSlack_[iColumn]+primalStep*change;2591 double gapProduct=dualValue*primalValue;2662 CoinWorkDouble change = rhsU_[iColumn]-deltaX_[iColumn]; 2663 CoinWorkDouble dualValue=wVec_[iColumn]+dualStep*deltaW_[iColumn]; 2664 CoinWorkDouble primalValue=upperSlack_[iColumn]+primalStep*change; 2665 CoinWorkDouble gapProduct=dualValue*primalValue; 2592 2666 if (gapProduct>0.0&&dualValue<0.0) 2593 2667 gapProduct = - gapProduct; … … 2598 2672 iColumn,primalValue,dualValue,gapProduct); 2599 2673 #endif 2600 double value= 0.0;2674 CoinWorkDouble value= 0.0; 2601 2675 if (gapProduct<minBeta) { 2602 2676 value= (minBeta-gapProduct); … … 2615 2689 if (cholesky_->type()<20) { 2616 2690 for (iColumn=0;iColumn<numberTotal;iColumn++) { 2617 double value = rhsC_[iColumn];2618 double zValue = rhsZ_[iColumn];2619 double wValue = rhsW_[iColumn];2691 CoinWorkDouble value = rhsC_[iColumn]; 2692 CoinWorkDouble zValue = rhsZ_[iColumn]; 2693 CoinWorkDouble wValue = rhsW_[iColumn]; 2620 2694 #if 0 2621 2695 #if 1 … … 2644 2718 #else 2645 2719 if (lowerBound(iColumn)) { 2646 double gHat = zValue + zVec_[iColumn]*rhsL_[iColumn];2720 CoinWorkDouble gHat = zValue + zVec_[iColumn]*rhsL_[iColumn]; 2647 2721 value -= gHat/(lowerSlack_[iColumn]+extra); 2648 2722 } 2649 2723 if (upperBound(iColumn)) { 2650 double hHat = wValue - wVec_[iColumn]*rhsU_[iColumn];2724 CoinWorkDouble hHat = wValue - wVec_[iColumn]*rhsU_[iColumn]; 2651 2725 value += hHat/(upperSlack_[iColumn]+extra); 2652 2726 } … … 2667 2741 // KKT 2668 2742 for (iColumn=0;iColumn<numberTotal;iColumn++) { 2669 double value = rhsC_[iColumn];2670 double zValue = rhsZ_[iColumn];2671 double wValue = rhsW_[iColumn];2743 CoinWorkDouble value = rhsC_[iColumn]; 2744 CoinWorkDouble zValue = rhsZ_[iColumn]; 2745 CoinWorkDouble wValue = rhsW_[iColumn]; 2672 2746 if (lowerBound(iColumn)) { 2673 double gHat = zValue + zVec_[iColumn]*rhsL_[iColumn];2747 CoinWorkDouble gHat = zValue + zVec_[iColumn]*rhsL_[iColumn]; 2674 2748 value -= gHat/(lowerSlack_[iColumn]+extra); 2675 2749 } 2676 2750 if (upperBound(iColumn)) { 2677 double hHat = wValue - wVec_[iColumn]*rhsU_[iColumn];2751 CoinWorkDouble hHat = wValue - wVec_[iColumn]*rhsU_[iColumn]; 2678 2752 value += hHat/(upperSlack_[iColumn]+extra); 2679 2753 } … … 2684 2758 //method: sees if looks plausible change in complementarity 2685 2759 bool ClpPredictorCorrector::checkGoodMove(const bool doCorrector, 2686 double & bestNextGap,2760 CoinWorkDouble & bestNextGap, 2687 2761 bool allowIncreasingGap) 2688 2762 { 2689 const double beta3 = 0.99997;2763 const CoinWorkDouble beta3 = 0.99997; 2690 2764 bool goodMove=false; 2691 2765 int nextNumber; 2692 2766 int nextNumberItems; 2693 2767 int numberTotal = numberRows_+numberColumns_; 2694 double returnGap=bestNextGap;2695 double nextGap=complementarityGap(nextNumber,nextNumberItems,2);2768 CoinWorkDouble returnGap=bestNextGap; 2769 CoinWorkDouble nextGap=complementarityGap(nextNumber,nextNumberItems,2); 2696 2770 #ifndef NO_RTTI 2697 2771 ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(objective_)); … … 2711 2785 returnGap=nextGap; 2712 2786 } 2713 double step;2787 CoinWorkDouble step; 2714 2788 if (actualDualStep_>actualPrimalStep_) { 2715 2789 step=actualDualStep_; … … 2717 2791 step=actualPrimalStep_; 2718 2792 } 2719 double testValue=1.0-step*(1.0-beta3);2793 CoinWorkDouble testValue=1.0-step*(1.0-beta3); 2720 2794 //testValue=0.0; 2721 2795 testValue*=complementarityGap_; … … 2729 2803 //step=actualPrimalStep_; 2730 2804 //} 2731 double gap = bestNextGap;2805 CoinWorkDouble gap = bestNextGap; 2732 2806 goodMove=checkGoodMove2(step,gap,allowIncreasingGap); 2733 2807 if (goodMove) … … 2760 2834 while (!goodMove) { 2761 2835 pass++; 2762 double gap = bestNextGap;2836 CoinWorkDouble gap = bestNextGap; 2763 2837 goodMove=checkGoodMove2(step,gap,allowIncreasingGap); 2764 2838 if (goodMove||pass>3) { … … 2793 2867 if (goodMove) { 2794 2868 //compute delta in objectives 2795 double deltaObjectivePrimal=0.0;2796 double deltaObjectiveDual=2869 CoinWorkDouble deltaObjectivePrimal=0.0; 2870 CoinWorkDouble deltaObjectiveDual= 2797 2871 innerProduct(deltaY_,numberRows_, 2798 2872 rhsFixRegion_); 2799 double error=0.0;2800 double * workArray = workArray_;2873 CoinWorkDouble error=0.0; 2874 CoinWorkDouble * workArray = workArray_; 2801 2875 CoinZeroN(workArray,numberColumns_); 2802 2876 CoinMemcpyN(deltaY_,numberRows_,workArray+numberColumns_); 2803 2877 matrix_->transposeTimes(-1.0,deltaY_,workArray); 2804 // double sumPerturbCost=0.0;2878 //CoinWorkDouble sumPerturbCost=0.0; 2805 2879 for (int iColumn=0;iColumn<numberTotal;iColumn++) { 2806 2880 if (!flagged(iColumn)) { … … 2813 2887 deltaObjectiveDual-=deltaW_[iColumn]*upper_[iColumn]; 2814 2888 } 2815 double change = fabs(workArray_[iColumn]-deltaZ_[iColumn]+deltaW_[iColumn]);2889 CoinWorkDouble change = CoinAbs(workArray_[iColumn]-deltaZ_[iColumn]+deltaW_[iColumn]); 2816 2890 error = CoinMax(change,error); 2817 2891 } … … 2819 2893 } 2820 2894 //deltaObjectivePrimal+=sumPerturbCost*linearPerturbation_; 2821 double testValue;2895 CoinWorkDouble testValue; 2822 2896 if (error>0.0) { 2823 2897 testValue=1.0e1*CoinMax(maximumDualError_,1.0e-12)/error; … … 2828 2902 if (testValue<actualDualStep_&&!quadraticObj) { 2829 2903 handler_->message(CLP_BARRIER_REDUCING,messages_) 2830 <<"dual"<<actualDualStep_ 2831 << testValue 2904 <<"dual"<<static_cast<double>(actualDualStep_) 2905 << static_cast<double>(testValue) 2832 2906 <<CoinMessageEol; 2833 2907 actualDualStep_=testValue; … … 2838 2912 //check change in AX not too much 2839 2913 //??? could be dropped row going infeasible 2840 double ratio = 1.0e1*CoinMax(maximumRHSError_,1.0e-12)/maximumRHSChange_;2914 CoinWorkDouble ratio = 1.0e1*CoinMax(maximumRHSError_,1.0e-12)/maximumRHSChange_; 2841 2915 if (ratio<actualPrimalStep_) { 2842 2916 handler_->message(CLP_BARRIER_REDUCING,messages_) 2843 <<"primal"<<actualPrimalStep_ 2844 <<ratio 2917 <<"primal"<<static_cast<double>(actualPrimalStep_) 2918 <<static_cast<double>(ratio) 2845 2919 <<CoinMessageEol; 2846 2920 if (ratio>1.0e-6) { … … 2857 2931 } 2858 2932 //: checks for one step size 2859 bool ClpPredictorCorrector::checkGoodMove2( double move,2860 double & bestNextGap,2933 bool ClpPredictorCorrector::checkGoodMove2(CoinWorkDouble move, 2934 CoinWorkDouble & bestNextGap, 2861 2935 bool allowIncreasingGap) 2862 2936 { 2863 double complementarityMultiplier =1.0/numberComplementarityPairs_; 2864 const double gamma = 1.0e-8; 2865 const double gammap = 1.0e-8; 2866 double gammad = 1.0e-8; 2937 CoinWorkDouble complementarityMultiplier =1.0/numberComplementarityPairs_; 2938 #if CLP_LONG_CHOLESKY<2 2939 const CoinWorkDouble gamma = 1.0e-8; 2940 const CoinWorkDouble gammap = 1.0e-8; 2941 CoinWorkDouble gammad = 1.0e-8; 2942 #else 2943 const CoinWorkDouble gamma = 1.0e-8; 2944 const CoinWorkDouble gammap = 1.0e-8; 2945 CoinWorkDouble gammad = 1.0e-8; 2946 #endif 2867 2947 int nextNumber; 2868 2948 int nextNumberItems; 2869 double nextGap=complementarityGap(nextNumber,nextNumberItems,2);2949 CoinWorkDouble nextGap=complementarityGap(nextNumber,nextNumberItems,2); 2870 2950 if (nextGap>bestNextGap&&!allowIncreasingGap) 2871 2951 return false; 2872 double lowerBoundGap = gamma*nextGap*complementarityMultiplier;2952 CoinWorkDouble lowerBoundGap = gamma*nextGap*complementarityMultiplier; 2873 2953 bool goodMove=true; 2874 2954 int iColumn; … … 2876 2956 if (!flagged(iColumn)) { 2877 2957 if (lowerBound(iColumn)) { 2878 double part1=lowerSlack_[iColumn]+actualPrimalStep_*deltaSL_[iColumn];2879 double part2=zVec_[iColumn]+actualDualStep_*deltaZ_[iColumn];2958 CoinWorkDouble part1=lowerSlack_[iColumn]+actualPrimalStep_*deltaSL_[iColumn]; 2959 CoinWorkDouble part2=zVec_[iColumn]+actualDualStep_*deltaZ_[iColumn]; 2880 2960 if (part1*part2<lowerBoundGap) { 2881 2961 goodMove=false; … … 2884 2964 } 2885 2965 if (upperBound(iColumn)) { 2886 double part1=upperSlack_[iColumn]+actualPrimalStep_*deltaSU_[iColumn];2887 double part2=wVec_[iColumn]+actualDualStep_*deltaW_[iColumn];2966 CoinWorkDouble part1=upperSlack_[iColumn]+actualPrimalStep_*deltaSU_[iColumn]; 2967 CoinWorkDouble part2=wVec_[iColumn]+actualDualStep_*deltaW_[iColumn]; 2888 2968 if (part1*part2<lowerBoundGap) { 2889 2969 goodMove=false; … … 2893 2973 } 2894 2974 } 2895 double * nextDj=NULL;2896 double maximumDualError = maximumDualError_;2975 CoinWorkDouble * nextDj=NULL; 2976 CoinWorkDouble maximumDualError = maximumDualError_; 2897 2977 #ifndef NO_RTTI 2898 2978 ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(objective_)); … … 2902 2982 quadraticObj = (static_cast< ClpQuadraticObjective*>(objective_)); 2903 2983 #endif 2984 CoinWorkDouble * dualArray = reinterpret_cast<CoinWorkDouble *>(dual_); 2904 2985 if (quadraticObj) { 2905 2986 // change gammad 2906 2987 gammad=1.0e-4; 2907 double gamma2 = gamma_*gamma_;2908 nextDj = new double [numberColumns_];2909 double * nextSolution = new double [numberColumns_];2988 CoinWorkDouble gamma2 = gamma_*gamma_; 2989 nextDj = new CoinWorkDouble [numberColumns_]; 2990 CoinWorkDouble * nextSolution = new CoinWorkDouble [numberColumns_]; 2910 2991 // put next primal into nextSolution 2911 2992 for ( iColumn=0;iColumn<numberColumns_;iColumn++) { … … 2919 3000 // do reduced costs 2920 3001 CoinMemcpyN(cost_,numberColumns_,nextDj); 2921 matrix_->transposeTimes(-1.0,dual _,nextDj);3002 matrix_->transposeTimes(-1.0,dualArray,nextDj); 2922 3003 matrix_->transposeTimes(-actualDualStep_,deltaY_,nextDj); 2923 3004 quadraticDjs(nextDj,nextSolution,1.0); … … 2927 3008 for (int iColumn=0;iColumn<numberColumns_;iColumn++) { 2928 3009 if (!fixedOrFree(iColumn)) { 2929 double newZ=0.0;2930 double newW=0.0;3010 CoinWorkDouble newZ=0.0; 3011 CoinWorkDouble newW=0.0; 2931 3012 if (lowerBound(iColumn)) { 2932 3013 newZ=zVec_[iColumn]+actualDualStep_*deltaZ_[iColumn]; … … 2936 3017 } 2937 3018 if (columnQuadraticLength[iColumn]) { 2938 double gammaTerm = gamma2;3019 CoinWorkDouble gammaTerm = gamma2; 2939 3020 if (primalR_) 2940 3021 gammaTerm += primalR_[iColumn]; 2941 // double dualInfeasibility=3022 //CoinWorkDouble dualInfeasibility= 2942 3023 //dj_[iColumn]-zVec_[iColumn]+wVec_[iColumn] 2943 3024 //+gammaTerm*solution_[iColumn]; 2944 double newInfeasibility=3025 CoinWorkDouble newInfeasibility= 2945 3026 nextDj[iColumn]-newZ+newW 2946 3027 +gammaTerm*(solution_[iColumn]+actualPrimalStep_*deltaX_[iColumn]); 2947 3028 maximumDualError = CoinMax(maximumDualError,newInfeasibility); 2948 //if ( fabs(newInfeasibility)>CoinMax(2000.0*maximumDualError_,1.0e-2)) {3029 //if (CoinAbs(newInfeasibility)>CoinMax(2000.0*maximumDualError_,1.0e-2)) { 2949 3030 //if (dualInfeasibility*newInfeasibility<0.0) { 2950 3031 // printf("%d current %g next %g\n",iColumn,dualInfeasibility, … … 2962 3043 solutionNorm_=rhsNorm_; 2963 3044 } 2964 double errorCheck=maximumRHSError_/solutionNorm_;3045 CoinWorkDouble errorCheck=maximumRHSError_/solutionNorm_; 2965 3046 if (errorCheck<maximumBoundInfeasibility_) { 2966 3047 errorCheck=maximumBoundInfeasibility_; … … 2987 3068 // updateSolution. Updates solution at end of iteration 2988 3069 //returns number fixed 2989 int ClpPredictorCorrector::updateSolution( double nextGap)3070 int ClpPredictorCorrector::updateSolution(CoinWorkDouble nextGap) 2990 3071 { 3072 CoinWorkDouble * dualArray = reinterpret_cast<CoinWorkDouble *>(dual_); 2991 3073 int numberTotal = numberRows_+numberColumns_; 2992 3074 //update pi 2993 multiplyAdd(deltaY_,numberRows_,actualDualStep_,dual _,1.0);3075 multiplyAdd(deltaY_,numberRows_,actualDualStep_,dualArray,1.0); 2994 3076 CoinZeroN(errorRegion_,numberRows_); 2995 3077 CoinZeroN(rhsFixRegion_,numberRows_); 2996 double maximumRhsInfeasibility=0.0;2997 double maximumBoundInfeasibility=0.0;2998 double maximumDualError=1.0e-12;2999 double primalObjectiveValue=0.0;3000 double dualObjectiveValue=0.0;3001 double solutionNorm=1.0e-12;3078 CoinWorkDouble maximumRhsInfeasibility=0.0; 3079 CoinWorkDouble maximumBoundInfeasibility=0.0; 3080 CoinWorkDouble maximumDualError=1.0e-12; 3081 CoinWorkDouble primalObjectiveValue=0.0; 3082 CoinWorkDouble dualObjectiveValue=0.0; 3083 CoinWorkDouble solutionNorm=1.0e-12; 3002 3084 int numberKilled=0; 3003 double freeMultiplier=1.0e6;3004 double trueNorm =diagonalNorm_/diagonalScaleFactor_;3085 CoinWorkDouble freeMultiplier=1.0e6; 3086 CoinWorkDouble trueNorm =diagonalNorm_/diagonalScaleFactor_; 3005 3087 if (freeMultiplier<trueNorm) { 3006 3088 freeMultiplier=trueNorm; … … 3010 3092 } 3011 3093 freeMultiplier=0.5/freeMultiplier; 3012 double condition = fabs(cholesky_->choleskyCondition());3094 CoinWorkDouble condition = CoinAbs(cholesky_->choleskyCondition()); 3013 3095 bool caution; 3014 3096 if ((condition<1.0e10&&trueNorm<1.0e12)||numberIterations_<20) { … … 3017 3099 caution=true; 3018 3100 } 3019 double extra=eExtra;3020 const double largeFactor=1.0e2;3021 double largeGap=largeFactor*solutionNorm_;3101 CoinWorkDouble extra=eExtra; 3102 const CoinWorkDouble largeFactor=1.0e2; 3103 CoinWorkDouble largeGap=largeFactor*solutionNorm_; 3022 3104 if (largeGap<largeFactor) { 3023 3105 largeGap=largeFactor; 3024 3106 } 3025 double dualFake=0.0;3026 double dualTolerance = dblParam_[ClpDualTolerance];3107 CoinWorkDouble dualFake=0.0; 3108 CoinWorkDouble dualTolerance = dblParam_[ClpDualTolerance]; 3027 3109 dualTolerance=dualTolerance/scaleFactor_; 3028 3110 if (dualTolerance<1.0e-12) { 3029 3111 dualTolerance=1.0e-12; 3030 3112 } 3031 double offsetObjective=0.0; 3032 const double killTolerance=primalTolerance(); 3033 double qDiagonal; 3113 CoinWorkDouble offsetObjective=0.0; 3114 const CoinWorkDouble killTolerance=primalTolerance(); 3115 CoinWorkDouble qDiagonal; 3116 #if CLP_LONG_CHOLESKY<2 3034 3117 if (mu_<1.0) { 3035 3118 qDiagona