Changeset 1371


Ignore:
Timestamp:
Jun 12, 2009 12:29:04 PM (10 years ago)
Author:
forrest
Message:

changes to try and make faster

Location:
trunk/Clp/src
Files:
14 edited

Legend:

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

    r1370 r1371  
    18161816#endif
    18171817#ifdef COIN_HAS_CLP
     1818  parameters[numberParameters++]=
     1819    CbcOrClpParam("fact!orization","Which factorization to use",
     1820                  "normal",FACTORIZATION);
     1821  parameters[numberParameters-1].append("dense");
     1822  parameters[numberParameters-1].append("simple");
     1823  parameters[numberParameters-1].append("osl");
     1824  parameters[numberParameters-1].setLonghelp
     1825    (
     1826     "The default is to use the normal CoinFactorization, but \
     1827other choices are a dense one, osl's or one designed for small problems."
     1828     );
    18181829  parameters[numberParameters++]=
    18191830    CbcOrClpParam("fakeB!ound","All bounds <= this value - DEBUG",
  • trunk/Clp/src/CbcOrClpParam.hpp

    r1370 r1371  
    7676    PRIMALPIVOT,PRESOLVE,CRASH,BIASLU,PERTURBATION,MESSAGES,AUTOSCALE,
    7777    CHOLESKY,KKT,BARRIERSCALE,GAMMA,CROSSOVER,PFI,INTPRINT,VECTOR,
     78    FACTORIZATION,
    7879   
    7980    NODESTRATEGY = 251,BRANCHSTRATEGY,CUTSSTRATEGY,HEURISTICSTRATEGY,
  • trunk/Clp/src/ClpCholeskyBase.cpp

    r1370 r1371  
    36933693/* Uses factorization to solve. */
    36943694void
    3695 ClpCholeskyBase::solve (double * region)
     3695ClpCholeskyBase::solve (CoinWorkDouble * region)
    36963696{
    36973697  if (!whichDense_) {
     
    37123712    }
    37133713    // solve
    3714 #if CLP_LONG_CHOLESKY>0
    3715     dense_->solveLong(change);
    3716 #else
    37173714    dense_->solve(change);
    3718 #endif
    37193715    for (i=0;i<numberDense;i++) {
    37203716      const longDouble * a = denseColumn_+i*numberRows_;
     
    37323728*/
    37333729void
    3734 ClpCholeskyBase::solve(double * region, int type)
     3730ClpCholeskyBase::solve(CoinWorkDouble * region, int type)
    37353731{
    37363732#ifdef CLP_DEBUG
     
    37903786      int nDense = numberRows_-firstDense_;
    37913787      dense.reserveSpace(this,nDense);
    3792 #if CLP_LONG_CHOLESKY!=1
    3793       dense.solveLong(work+firstDense_);
    3794 #else
    3795       dense.solveLongWork(work+firstDense_);
    3796 #endif
     3788      dense.solve(work+firstDense_);
    37973789      for (i=numberRows_-1;i>=firstDense_;i--) {
    37983790        CoinWorkDouble value=work[i];
     
    38413833      int nDense = numberRows_-firstDense_;
    38423834      dense.reserveSpace(this,nDense);
    3843       dense.solveLong(work+firstDense_);
     3835      dense.solve(work+firstDense_);
    38443836      for (i=numberRows_-1;i>=firstDense_;i--) {
    38453837        CoinWorkDouble value=work[i];
     
    38713863#endif
    38723864}
    3873 #if CLP_LONG_CHOLESKY
     3865#if 0 //CLP_LONG_CHOLESKY
    38743866/* Uses factorization to solve. */
    38753867void
     
    38983890    int nDense = numberRows_-firstDense_;
    38993891    dense.reserveSpace(this,nDense);
    3900 #if CLP_LONG_CHOLESKY!=1
    3901     dense.solveLong(work+firstDense_);
    3902 #else
    3903     dense.solveLongWork(work+firstDense_);
    3904 #endif
     3892    dense.solve(work+firstDense_);
    39053893    for (i=numberRows_-1;i>=firstDense_;i--) {
    39063894      CoinWorkDouble value=work[i];
  • trunk/Clp/src/ClpCholeskyBase.hpp

    r1370 r1371  
    7070  virtual int factorize(const CoinWorkDouble * diagonal, int * rowsDropped) ;
    7171  /** Uses factorization to solve. */
    72   virtual void solve (double * region) ;
     72  virtual void solve (CoinWorkDouble * region) ;
    7373  /** Uses factorization to solve. - given as if KKT.
    7474   region1 is rows+columns, region2 is rows */
    7575  virtual void solveKKT (CoinWorkDouble * region1, CoinWorkDouble * region2, const CoinWorkDouble * diagonal,
    7676                         CoinWorkDouble diagonalScaleFactor);
    77 #if CLP_LONG_CHOLESKY
    78   /** Uses factorization to solve. */
    79   virtual void solve (CoinWorkDouble * region) ;
    80 #endif
    8177private:
    8278  /// AMD ordering
     
    195191  If 1 and 2 then diagonal has sqrt of inverse otherwise inverse
    196192  */
    197   void solve(double * region, int type);
     193  void solve(CoinWorkDouble * region, int type);
    198194  /// Forms ADAT - returns nonzero if not enough memory
    199195  int preOrder(bool lowerTriangular, bool includeDiagonal, bool doKKT);
  • trunk/Clp/src/ClpCholeskyDense.cpp

    r1370 r1371  
    150150/* Factorize - filling in rowsDropped and returning number dropped */
    151151int
    152 ClpCholeskyDense::factorize(const double * diagonal, int * rowsDropped)
     152ClpCholeskyDense::factorize(const CoinWorkDouble * diagonal, int * rowsDropped)
    153153{
    154154  assert (!borrowSpace_);
     
    184184    work--; /* skip diagonal*/
    185185    int addOffset=numberRows_-1;
    186     const double * diagonalSlack = diagonal + numberColumns;
     186    const CoinWorkDouble * diagonalSlack = diagonal + numberColumns;
    187187    /* largest in initial matrix*/
    188188    CoinWorkDouble largest2=1.0e-20;
     
    12151215/* Uses factorization to solve. */
    12161216void
    1217 ClpCholeskyDense::solve (double * region)
     1217ClpCholeskyDense::solve (CoinWorkDouble * region)
    12181218{
    12191219#ifdef CHOL_COMPARE
     
    13201320/* Forward part of solve 1*/
    13211321void
    1322 ClpCholeskyDense::solveF1(longDouble * a,int n,double * region)
     1322ClpCholeskyDense::solveF1(longDouble * a,int n,CoinWorkDouble * region)
    13231323{
    13241324  int j, k;
     
    13351335/* Forward part of solve 2*/
    13361336void
    1337 ClpCholeskyDense::solveF2(longDouble * a,int n,double * region, double * region2)
     1337ClpCholeskyDense::solveF2(longDouble * a,int n,CoinWorkDouble * region, CoinWorkDouble * region2)
    13381338{
    13391339  int j, k;
     
    14471447/* Backward part of solve 1*/
    14481448void
    1449 ClpCholeskyDense::solveB1(longDouble * a,int n,double * region)
     1449ClpCholeskyDense::solveB1(longDouble * a,int n,CoinWorkDouble * region)
    14501450{
    14511451  int j, k;
     
    14621462/* Backward part of solve 2*/
    14631463void
    1464 ClpCholeskyDense::solveB2(longDouble * a,int n,double * region, double * region2)
     1464ClpCholeskyDense::solveB2(longDouble * a,int n,CoinWorkDouble * region, CoinWorkDouble * region2)
    14651465{
    14661466  int j, k;
     
    15731573#endif
    15741574}
    1575 /* Uses factorization to solve. */
    1576 void
    1577 ClpCholeskyDense::solveLong (CoinWorkDouble * region)
    1578 {
    1579   int numberBlocks = (numberRows_+BLOCK-1)>>BLOCKSHIFT;
    1580   /* later align on boundary*/
    1581   longDouble * a = sparseFactor_+BLOCKSQ*numberBlocks;
    1582   int iBlock;
    1583   longDouble * aa = a;
    1584   int iColumn;
    1585   for (iBlock=0;iBlock<numberBlocks;iBlock++) {
    1586     int nChunk;
    1587     int jBlock;
    1588     int iDo = iBlock*BLOCK;
    1589     int base=iDo;
    1590     if (iDo+BLOCK>numberRows_) {
    1591       nChunk=numberRows_-iDo;
    1592     } else {
    1593       nChunk=BLOCK;
    1594     }
    1595     solveF1Long(aa,nChunk,region+iDo);
    1596     for (jBlock=iBlock+1;jBlock<numberBlocks;jBlock++) {
    1597       base+=BLOCK;
    1598       aa+=BLOCKSQ;
    1599       if (base+BLOCK>numberRows_) {
    1600         nChunk=numberRows_-base;
    1601       } else {
    1602         nChunk=BLOCK;
    1603       }
    1604       solveF2Long(aa,nChunk,region+iDo,region+base);
    1605     }
    1606     aa+=BLOCKSQ;
    1607   }
    1608   /* do diagonal outside*/
    1609   for (iColumn=0;iColumn<numberRows_;iColumn++)
    1610     region[iColumn] *= diagonal_[iColumn];
    1611   int offset=((numberBlocks*(numberBlocks+1))>>1);
    1612   aa = a+number_entries(offset-1);
    1613   int lBase=(numberBlocks-1)*BLOCK;
    1614   for (iBlock=numberBlocks-1;iBlock>=0;iBlock--) {
    1615     int nChunk;
    1616     int jBlock;
    1617     int triBase=iBlock*BLOCK;
    1618     int iBase=lBase;
    1619     for (jBlock=iBlock+1;jBlock<numberBlocks;jBlock++) {
    1620       if (iBase+BLOCK>numberRows_) {
    1621         nChunk=numberRows_-iBase;
    1622       } else {
    1623         nChunk=BLOCK;
    1624       }
    1625       solveB2Long(aa,nChunk,region+triBase,region+iBase);
    1626       iBase-=BLOCK;
    1627       aa-=BLOCKSQ;
    1628     }
    1629     if (triBase+BLOCK>numberRows_) {
    1630       nChunk=numberRows_-triBase;
    1631     } else {
    1632       nChunk=BLOCK;
    1633     }
    1634     solveB1Long(aa,nChunk,region+triBase);
    1635     aa-=BLOCKSQ;
    1636   }
    1637 }
    1638 /* Forward part of solve 1*/
    1639 void
    1640 ClpCholeskyDense::solveF1Long(longDouble * a,int n,CoinWorkDouble * region)
    1641 {
    1642   int j, k;
    1643   CoinWorkDouble t00;
    1644   for (j = 0; j < n; j ++) {
    1645     t00 = region[j];
    1646     for (k = 0; k < j; ++k) {
    1647       t00 -= region[k] * a[j + k * BLOCK];
    1648     }
    1649     /*t00*=a[j + j * BLOCK];*/
    1650     region[j] = t00;
    1651   }
    1652 }
    1653 /* Forward part of solve 2*/
    1654 void
    1655 ClpCholeskyDense::solveF2Long(longDouble * a,int n,CoinWorkDouble * region, CoinWorkDouble * region2)
    1656 {
    1657   int j, k;
    1658 #ifdef BLOCKUNROLL
    1659   if (n==BLOCK) {
    1660     for (k = 0; k < BLOCK; k+=4) {
    1661       CoinWorkDouble t0 = region2[0];
    1662       CoinWorkDouble t1 = region2[1];
    1663       CoinWorkDouble t2 = region2[2];
    1664       CoinWorkDouble t3 = region2[3];
    1665       t0 -= region[0] * a[0 + 0 * BLOCK];
    1666       t1 -= region[0] * a[1 + 0 * BLOCK];
    1667       t2 -= region[0] * a[2 + 0 * BLOCK];
    1668       t3 -= region[0] * a[3 + 0 * BLOCK];
    1669 
    1670       t0 -= region[1] * a[0 + 1 * BLOCK];
    1671       t1 -= region[1] * a[1 + 1 * BLOCK];
    1672       t2 -= region[1] * a[2 + 1 * BLOCK];
    1673       t3 -= region[1] * a[3 + 1 * BLOCK];
    1674 
    1675       t0 -= region[2] * a[0 + 2 * BLOCK];
    1676       t1 -= region[2] * a[1 + 2 * BLOCK];
    1677       t2 -= region[2] * a[2 + 2 * BLOCK];
    1678       t3 -= region[2] * a[3 + 2 * BLOCK];
    1679 
    1680       t0 -= region[3] * a[0 + 3 * BLOCK];
    1681       t1 -= region[3] * a[1 + 3 * BLOCK];
    1682       t2 -= region[3] * a[2 + 3 * BLOCK];
    1683       t3 -= region[3] * a[3 + 3 * BLOCK];
    1684 
    1685       t0 -= region[4] * a[0 + 4 * BLOCK];
    1686       t1 -= region[4] * a[1 + 4 * BLOCK];
    1687       t2 -= region[4] * a[2 + 4 * BLOCK];
    1688       t3 -= region[4] * a[3 + 4 * BLOCK];
    1689 
    1690       t0 -= region[5] * a[0 + 5 * BLOCK];
    1691       t1 -= region[5] * a[1 + 5 * BLOCK];
    1692       t2 -= region[5] * a[2 + 5 * BLOCK];
    1693       t3 -= region[5] * a[3 + 5 * BLOCK];
    1694 
    1695       t0 -= region[6] * a[0 + 6 * BLOCK];
    1696       t1 -= region[6] * a[1 + 6 * BLOCK];
    1697       t2 -= region[6] * a[2 + 6 * BLOCK];
    1698       t3 -= region[6] * a[3 + 6 * BLOCK];
    1699 
    1700       t0 -= region[7] * a[0 + 7 * BLOCK];
    1701       t1 -= region[7] * a[1 + 7 * BLOCK];
    1702       t2 -= region[7] * a[2 + 7 * BLOCK];
    1703       t3 -= region[7] * a[3 + 7 * BLOCK];
    1704 #if BLOCK>8
    1705       t0 -= region[8] * a[0 + 8 * BLOCK];
    1706       t1 -= region[8] * a[1 + 8 * BLOCK];
    1707       t2 -= region[8] * a[2 + 8 * BLOCK];
    1708       t3 -= region[8] * a[3 + 8 * BLOCK];
    1709 
    1710       t0 -= region[9] * a[0 + 9 * BLOCK];
    1711       t1 -= region[9] * a[1 + 9 * BLOCK];
    1712       t2 -= region[9] * a[2 + 9 * BLOCK];
    1713       t3 -= region[9] * a[3 + 9 * BLOCK];
    1714 
    1715       t0 -= region[10] * a[0 + 10 * BLOCK];
    1716       t1 -= region[10] * a[1 + 10 * BLOCK];
    1717       t2 -= region[10] * a[2 + 10 * BLOCK];
    1718       t3 -= region[10] * a[3 + 10 * BLOCK];
    1719 
    1720       t0 -= region[11] * a[0 + 11 * BLOCK];
    1721       t1 -= region[11] * a[1 + 11 * BLOCK];
    1722       t2 -= region[11] * a[2 + 11 * BLOCK];
    1723       t3 -= region[11] * a[3 + 11 * BLOCK];
    1724 
    1725       t0 -= region[12] * a[0 + 12 * BLOCK];
    1726       t1 -= region[12] * a[1 + 12 * BLOCK];
    1727       t2 -= region[12] * a[2 + 12 * BLOCK];
    1728       t3 -= region[12] * a[3 + 12 * BLOCK];
    1729 
    1730       t0 -= region[13] * a[0 + 13 * BLOCK];
    1731       t1 -= region[13] * a[1 + 13 * BLOCK];
    1732       t2 -= region[13] * a[2 + 13 * BLOCK];
    1733       t3 -= region[13] * a[3 + 13 * BLOCK];
    1734 
    1735       t0 -= region[14] * a[0 + 14 * BLOCK];
    1736       t1 -= region[14] * a[1 + 14 * BLOCK];
    1737       t2 -= region[14] * a[2 + 14 * BLOCK];
    1738       t3 -= region[14] * a[3 + 14 * BLOCK];
    1739 
    1740       t0 -= region[15] * a[0 + 15 * BLOCK];
    1741       t1 -= region[15] * a[1 + 15 * BLOCK];
    1742       t2 -= region[15] * a[2 + 15 * BLOCK];
    1743       t3 -= region[15] * a[3 + 15 * BLOCK];
    1744 #endif
    1745       region2[0] = t0;
    1746       region2[1] = t1;
    1747       region2[2] = t2;
    1748       region2[3] = t3;
    1749       region2+=4;
    1750       a+=4;
    1751     }
    1752   } else {
    1753 #endif
    1754     for (k = 0; k < n; ++k) {
    1755       CoinWorkDouble t00 = region2[k];
    1756       for (j = 0; j < BLOCK; j ++) {
    1757         t00 -= region[j] * a[k + j * BLOCK];
    1758       }
    1759       region2[k] = t00;
    1760     }
    1761 #ifdef BLOCKUNROLL
    1762   }
    1763 #endif
    1764 }
    1765 /* Backward part of solve 1*/
    1766 void
    1767 ClpCholeskyDense::solveB1Long(longDouble * a,int n,CoinWorkDouble * region)
    1768 {
    1769   int j, k;
    1770   CoinWorkDouble t00;
    1771   for (j = n-1; j >=0; j --) {
    1772     t00 = region[j];
    1773     for (k = j+1; k < n; ++k) {
    1774       t00 -= region[k] * a[k + j * BLOCK];
    1775     }
    1776     /*t00*=a[j + j * BLOCK];*/
    1777     region[j] = t00;
    1778   }
    1779 }
    1780 /* Backward part of solve 2*/
    1781 void
    1782 ClpCholeskyDense::solveB2Long(longDouble * a,int n,CoinWorkDouble * region, CoinWorkDouble * region2)
    1783 {
    1784   int j, k;
    1785 #ifdef BLOCKUNROLL
    1786   if (n==BLOCK) {
    1787     for (j = 0; j < BLOCK; j +=4) {
    1788       CoinWorkDouble t0 = region[0];
    1789       CoinWorkDouble t1 = region[1];
    1790       CoinWorkDouble t2 = region[2];
    1791       CoinWorkDouble t3 = region[3];
    1792       t0 -= region2[0] * a[0 + 0*BLOCK];
    1793       t1 -= region2[0] * a[0 + 1*BLOCK];
    1794       t2 -= region2[0] * a[0 + 2*BLOCK];
    1795       t3 -= region2[0] * a[0 + 3*BLOCK];
    1796 
    1797       t0 -= region2[1] * a[1 + 0*BLOCK];
    1798       t1 -= region2[1] * a[1 + 1*BLOCK];
    1799       t2 -= region2[1] * a[1 + 2*BLOCK];
    1800       t3 -= region2[1] * a[1 + 3*BLOCK];
    1801 
    1802       t0 -= region2[2] * a[2 + 0*BLOCK];
    1803       t1 -= region2[2] * a[2 + 1*BLOCK];
    1804       t2 -= region2[2] * a[2 + 2*BLOCK];
    1805       t3 -= region2[2] * a[2 + 3*BLOCK];
    1806 
    1807       t0 -= region2[3] * a[3 + 0*BLOCK];
    1808       t1 -= region2[3] * a[3 + 1*BLOCK];
    1809       t2 -= region2[3] * a[3 + 2*BLOCK];
    1810       t3 -= region2[3] * a[3 + 3*BLOCK];
    1811 
    1812       t0 -= region2[4] * a[4 + 0*BLOCK];
    1813       t1 -= region2[4] * a[4 + 1*BLOCK];
    1814       t2 -= region2[4] * a[4 + 2*BLOCK];
    1815       t3 -= region2[4] * a[4 + 3*BLOCK];
    1816 
    1817       t0 -= region2[5] * a[5 + 0*BLOCK];
    1818       t1 -= region2[5] * a[5 + 1*BLOCK];
    1819       t2 -= region2[5] * a[5 + 2*BLOCK];
    1820       t3 -= region2[5] * a[5 + 3*BLOCK];
    1821 
    1822       t0 -= region2[6] * a[6 + 0*BLOCK];
    1823       t1 -= region2[6] * a[6 + 1*BLOCK];
    1824       t2 -= region2[6] * a[6 + 2*BLOCK];
    1825       t3 -= region2[6] * a[6 + 3*BLOCK];
    1826 
    1827       t0 -= region2[7] * a[7 + 0*BLOCK];
    1828       t1 -= region2[7] * a[7 + 1*BLOCK];
    1829       t2 -= region2[7] * a[7 + 2*BLOCK];
    1830       t3 -= region2[7] * a[7 + 3*BLOCK];
    1831 #if BLOCK>8
    1832 
    1833       t0 -= region2[8] * a[8 + 0*BLOCK];
    1834       t1 -= region2[8] * a[8 + 1*BLOCK];
    1835       t2 -= region2[8] * a[8 + 2*BLOCK];
    1836       t3 -= region2[8] * a[8 + 3*BLOCK];
    1837 
    1838       t0 -= region2[9] * a[9 + 0*BLOCK];
    1839       t1 -= region2[9] * a[9 + 1*BLOCK];
    1840       t2 -= region2[9] * a[9 + 2*BLOCK];
    1841       t3 -= region2[9] * a[9 + 3*BLOCK];
    1842 
    1843       t0 -= region2[10] * a[10 + 0*BLOCK];
    1844       t1 -= region2[10] * a[10 + 1*BLOCK];
    1845       t2 -= region2[10] * a[10 + 2*BLOCK];
    1846       t3 -= region2[10] * a[10 + 3*BLOCK];
    1847 
    1848       t0 -= region2[11] * a[11 + 0*BLOCK];
    1849       t1 -= region2[11] * a[11 + 1*BLOCK];
    1850       t2 -= region2[11] * a[11 + 2*BLOCK];
    1851       t3 -= region2[11] * a[11 + 3*BLOCK];
    1852 
    1853       t0 -= region2[12] * a[12 + 0*BLOCK];
    1854       t1 -= region2[12] * a[12 + 1*BLOCK];
    1855       t2 -= region2[12] * a[12 + 2*BLOCK];
    1856       t3 -= region2[12] * a[12 + 3*BLOCK];
    1857 
    1858       t0 -= region2[13] * a[13 + 0*BLOCK];
    1859       t1 -= region2[13] * a[13 + 1*BLOCK];
    1860       t2 -= region2[13] * a[13 + 2*BLOCK];
    1861       t3 -= region2[13] * a[13 + 3*BLOCK];
    1862 
    1863       t0 -= region2[14] * a[14 + 0*BLOCK];
    1864       t1 -= region2[14] * a[14 + 1*BLOCK];
    1865       t2 -= region2[14] * a[14 + 2*BLOCK];
    1866       t3 -= region2[14] * a[14 + 3*BLOCK];
    1867 
    1868       t0 -= region2[15] * a[15 + 0*BLOCK];
    1869       t1 -= region2[15] * a[15 + 1*BLOCK];
    1870       t2 -= region2[15] * a[15 + 2*BLOCK];
    1871       t3 -= region2[15] * a[15 + 3*BLOCK];
    1872 #endif
    1873       region[0] = t0;
    1874       region[1] = t1;
    1875       region[2] = t2;
    1876       region[3] = t3;
    1877       a+=4*BLOCK;
    1878       region+=4;
    1879     }
    1880   } else {
    1881 #endif
    1882     for (j = 0; j < BLOCK; j ++) {
    1883       CoinWorkDouble t00 = region[j];
    1884       for (k = 0; k < n; ++k) {
    1885         t00 -= region2[k] * a[k + j * BLOCK];
    1886       }
    1887       region[j] = t00;
    1888     }
    1889 #ifdef BLOCKUNROLL
    1890   }
    1891 #endif
    1892 }
    1893 /* Uses factorization to solve. */
    1894 void
    1895 ClpCholeskyDense::solveLongWork (CoinWorkDouble * region)
    1896 {
    1897   int numberBlocks = (numberRows_+BLOCK-1)>>BLOCKSHIFT;
    1898   /* later align on boundary*/
    1899   longDouble * a = sparseFactor_+BLOCKSQ*numberBlocks;
    1900   int iBlock;
    1901   longDouble * aa = a;
    1902   int iColumn;
    1903   for (iBlock=0;iBlock<numberBlocks;iBlock++) {
    1904     int nChunk;
    1905     int jBlock;
    1906     int iDo = iBlock*BLOCK;
    1907     int base=iDo;
    1908     if (iDo+BLOCK>numberRows_) {
    1909       nChunk=numberRows_-iDo;
    1910     } else {
    1911       nChunk=BLOCK;
    1912     }
    1913     solveF1LongWork(aa,nChunk,region+iDo);
    1914     for (jBlock=iBlock+1;jBlock<numberBlocks;jBlock++) {
    1915       base+=BLOCK;
    1916       aa+=BLOCKSQ;
    1917       if (base+BLOCK>numberRows_) {
    1918         nChunk=numberRows_-base;
    1919       } else {
    1920         nChunk=BLOCK;
    1921       }
    1922       solveF2LongWork(aa,nChunk,region+iDo,region+base);
    1923     }
    1924     aa+=BLOCKSQ;
    1925   }
    1926   /* do diagonal outside*/
    1927   for (iColumn=0;iColumn<numberRows_;iColumn++)
    1928     region[iColumn] *= diagonal_[iColumn];
    1929   int offset=((numberBlocks*(numberBlocks+1))>>1);
    1930   aa = a+number_entries(offset-1);
    1931   int lBase=(numberBlocks-1)*BLOCK;
    1932   for (iBlock=numberBlocks-1;iBlock>=0;iBlock--) {
    1933     int nChunk;
    1934     int jBlock;
    1935     int triBase=iBlock*BLOCK;
    1936     int iBase=lBase;
    1937     for (jBlock=iBlock+1;jBlock<numberBlocks;jBlock++) {
    1938       if (iBase+BLOCK>numberRows_) {
    1939         nChunk=numberRows_-iBase;
    1940       } else {
    1941         nChunk=BLOCK;
    1942       }
    1943       solveB2LongWork(aa,nChunk,region+triBase,region+iBase);
    1944       iBase-=BLOCK;
    1945       aa-=BLOCKSQ;
    1946     }
    1947     if (triBase+BLOCK>numberRows_) {
    1948       nChunk=numberRows_-triBase;
    1949     } else {
    1950       nChunk=BLOCK;
    1951     }
    1952     solveB1LongWork(aa,nChunk,region+triBase);
    1953     aa-=BLOCKSQ;
    1954   }
    1955 }
    1956 /* Forward part of solve 1*/
    1957 void
    1958 ClpCholeskyDense::solveF1LongWork(longDouble * a,int n,CoinWorkDouble * region)
    1959 {
    1960   int j, k;
    1961   CoinWorkDouble t00;
    1962   for (j = 0; j < n; j ++) {
    1963     t00 = region[j];
    1964     for (k = 0; k < j; ++k) {
    1965       t00 -= region[k] * a[j + k * BLOCK];
    1966     }
    1967     /*t00*=a[j + j * BLOCK];*/
    1968     region[j] = t00;
    1969   }
    1970 }
    1971 /* Forward part of solve 2*/
    1972 void
    1973 ClpCholeskyDense::solveF2LongWork(longDouble * a,int n,CoinWorkDouble * region, CoinWorkDouble * region2)
    1974 {
    1975   int j, k;
    1976 #ifdef BLOCKUNROLL
    1977   if (n==BLOCK) {
    1978     for (k = 0; k < BLOCK; k+=4) {
    1979       CoinWorkDouble t0 = region2[0];
    1980       CoinWorkDouble t1 = region2[1];
    1981       CoinWorkDouble t2 = region2[2];
    1982       CoinWorkDouble t3 = region2[3];
    1983       t0 -= region[0] * a[0 + 0 * BLOCK];
    1984       t1 -= region[0] * a[1 + 0 * BLOCK];
    1985       t2 -= region[0] * a[2 + 0 * BLOCK];
    1986       t3 -= region[0] * a[3 + 0 * BLOCK];
    1987 
    1988       t0 -= region[1] * a[0 + 1 * BLOCK];
    1989       t1 -= region[1] * a[1 + 1 * BLOCK];
    1990       t2 -= region[1] * a[2 + 1 * BLOCK];
    1991       t3 -= region[1] * a[3 + 1 * BLOCK];
    1992 
    1993       t0 -= region[2] * a[0 + 2 * BLOCK];
    1994       t1 -= region[2] * a[1 + 2 * BLOCK];
    1995       t2 -= region[2] * a[2 + 2 * BLOCK];
    1996       t3 -= region[2] * a[3 + 2 * BLOCK];
    1997 
    1998       t0 -= region[3] * a[0 + 3 * BLOCK];
    1999       t1 -= region[3] * a[1 + 3 * BLOCK];
    2000       t2 -= region[3] * a[2 + 3 * BLOCK];
    2001       t3 -= region[3] * a[3 + 3 * BLOCK];
    2002 
    2003       t0 -= region[4] * a[0 + 4 * BLOCK];
    2004       t1 -= region[4] * a[1 + 4 * BLOCK];
    2005       t2 -= region[4] * a[2 + 4 * BLOCK];
    2006       t3 -= region[4] * a[3 + 4 * BLOCK];
    2007 
    2008       t0 -= region[5] * a[0 + 5 * BLOCK];
    2009       t1 -= region[5] * a[1 + 5 * BLOCK];
    2010       t2 -= region[5] * a[2 + 5 * BLOCK];
    2011       t3 -= region[5] * a[3 + 5 * BLOCK];
    2012 
    2013       t0 -= region[6] * a[0 + 6 * BLOCK];
    2014       t1 -= region[6] * a[1 + 6 * BLOCK];
    2015       t2 -= region[6] * a[2 + 6 * BLOCK];
    2016       t3 -= region[6] * a[3 + 6 * BLOCK];
    2017 
    2018       t0 -= region[7] * a[0 + 7 * BLOCK];
    2019       t1 -= region[7] * a[1 + 7 * BLOCK];
    2020       t2 -= region[7] * a[2 + 7 * BLOCK];
    2021       t3 -= region[7] * a[3 + 7 * BLOCK];
    2022 #if BLOCK>8
    2023       t0 -= region[8] * a[0 + 8 * BLOCK];
    2024       t1 -= region[8] * a[1 + 8 * BLOCK];
    2025       t2 -= region[8] * a[2 + 8 * BLOCK];
    2026       t3 -= region[8] * a[3 + 8 * BLOCK];
    2027 
    2028       t0 -= region[9] * a[0 + 9 * BLOCK];
    2029       t1 -= region[9] * a[1 + 9 * BLOCK];
    2030       t2 -= region[9] * a[2 + 9 * BLOCK];
    2031       t3 -= region[9] * a[3 + 9 * BLOCK];
    2032 
    2033       t0 -= region[10] * a[0 + 10 * BLOCK];
    2034       t1 -= region[10] * a[1 + 10 * BLOCK];
    2035       t2 -= region[10] * a[2 + 10 * BLOCK];
    2036       t3 -= region[10] * a[3 + 10 * BLOCK];
    2037 
    2038       t0 -= region[11] * a[0 + 11 * BLOCK];
    2039       t1 -= region[11] * a[1 + 11 * BLOCK];
    2040       t2 -= region[11] * a[2 + 11 * BLOCK];
    2041       t3 -= region[11] * a[3 + 11 * BLOCK];
    2042 
    2043       t0 -= region[12] * a[0 + 12 * BLOCK];
    2044       t1 -= region[12] * a[1 + 12 * BLOCK];
    2045       t2 -= region[12] * a[2 + 12 * BLOCK];
    2046       t3 -= region[12] * a[3 + 12 * BLOCK];
    2047 
    2048       t0 -= region[13] * a[0 + 13 * BLOCK];
    2049       t1 -= region[13] * a[1 + 13 * BLOCK];
    2050       t2 -= region[13] * a[2 + 13 * BLOCK];
    2051       t3 -= region[13] * a[3 + 13 * BLOCK];
    2052 
    2053       t0 -= region[14] * a[0 + 14 * BLOCK];
    2054       t1 -= region[14] * a[1 + 14 * BLOCK];
    2055       t2 -= region[14] * a[2 + 14 * BLOCK];
    2056       t3 -= region[14] * a[3 + 14 * BLOCK];
    2057 
    2058       t0 -= region[15] * a[0 + 15 * BLOCK];
    2059       t1 -= region[15] * a[1 + 15 * BLOCK];
    2060       t2 -= region[15] * a[2 + 15 * BLOCK];
    2061       t3 -= region[15] * a[3 + 15 * BLOCK];
    2062 #endif
    2063       region2[0] = t0;
    2064       region2[1] = t1;
    2065       region2[2] = t2;
    2066       region2[3] = t3;
    2067       region2+=4;
    2068       a+=4;
    2069     }
    2070   } else {
    2071 #endif
    2072     for (k = 0; k < n; ++k) {
    2073       CoinWorkDouble t00 = region2[k];
    2074       for (j = 0; j < BLOCK; j ++) {
    2075         t00 -= region[j] * a[k + j * BLOCK];
    2076       }
    2077       region2[k] = t00;
    2078     }
    2079 #ifdef BLOCKUNROLL
    2080   }
    2081 #endif
    2082 }
    2083 /* Backward part of solve 1*/
    2084 void
    2085 ClpCholeskyDense::solveB1LongWork(longDouble * a,int n,CoinWorkDouble * region)
    2086 {
    2087   int j, k;
    2088   CoinWorkDouble t00;
    2089   for (j = n-1; j >=0; j --) {
    2090     t00 = region[j];
    2091     for (k = j+1; k < n; ++k) {
    2092       t00 -= region[k] * a[k + j * BLOCK];
    2093     }
    2094     /*t00*=a[j + j * BLOCK];*/
    2095     region[j] = t00;
    2096   }
    2097 }
    2098 /* Backward part of solve 2*/
    2099 void
    2100 ClpCholeskyDense::solveB2LongWork(longDouble * a,int n,CoinWorkDouble * region, CoinWorkDouble * region2)
    2101 {
    2102   int j, k;
    2103 #ifdef BLOCKUNROLL
    2104   if (n==BLOCK) {
    2105     for (j = 0; j < BLOCK; j +=4) {
    2106       CoinWorkDouble t0 = region[0];
    2107       CoinWorkDouble t1 = region[1];
    2108       CoinWorkDouble t2 = region[2];
    2109       CoinWorkDouble t3 = region[3];
    2110       t0 -= region2[0] * a[0 + 0*BLOCK];
    2111       t1 -= region2[0] * a[0 + 1*BLOCK];
    2112       t2 -= region2[0] * a[0 + 2*BLOCK];
    2113       t3 -= region2[0] * a[0 + 3*BLOCK];
    2114 
    2115       t0 -= region2[1] * a[1 + 0*BLOCK];
    2116       t1 -= region2[1] * a[1 + 1*BLOCK];
    2117       t2 -= region2[1] * a[1 + 2*BLOCK];
    2118       t3 -= region2[1] * a[1 + 3*BLOCK];
    2119 
    2120       t0 -= region2[2] * a[2 + 0*BLOCK];
    2121       t1 -= region2[2] * a[2 + 1*BLOCK];
    2122       t2 -= region2[2] * a[2 + 2*BLOCK];
    2123       t3 -= region2[2] * a[2 + 3*BLOCK];
    2124 
    2125       t0 -= region2[3] * a[3 + 0*BLOCK];
    2126       t1 -= region2[3] * a[3 + 1*BLOCK];
    2127       t2 -= region2[3] * a[3 + 2*BLOCK];
    2128       t3 -= region2[3] * a[3 + 3*BLOCK];
    2129 
    2130       t0 -= region2[4] * a[4 + 0*BLOCK];
    2131       t1 -= region2[4] * a[4 + 1*BLOCK];
    2132       t2 -= region2[4] * a[4 + 2*BLOCK];
    2133       t3 -= region2[4] * a[4 + 3*BLOCK];
    2134 
    2135       t0 -= region2[5] * a[5 + 0*BLOCK];
    2136       t1 -= region2[5] * a[5 + 1*BLOCK];
    2137       t2 -= region2[5] * a[5 + 2*BLOCK];
    2138       t3 -= region2[5] * a[5 + 3*BLOCK];
    2139 
    2140       t0 -= region2[6] * a[6 + 0*BLOCK];
    2141       t1 -= region2[6] * a[6 + 1*BLOCK];
    2142       t2 -= region2[6] * a[6 + 2*BLOCK];
    2143       t3 -= region2[6] * a[6 + 3*BLOCK];
    2144 
    2145       t0 -= region2[7] * a[7 + 0*BLOCK];
    2146       t1 -= region2[7] * a[7 + 1*BLOCK];
    2147       t2 -= region2[7] * a[7 + 2*BLOCK];
    2148       t3 -= region2[7] * a[7 + 3*BLOCK];
    2149 #if BLOCK>8
    2150 
    2151       t0 -= region2[8] * a[8 + 0*BLOCK];
    2152       t1 -= region2[8] * a[8 + 1*BLOCK];
    2153       t2 -= region2[8] * a[8 + 2*BLOCK];
    2154       t3 -= region2[8] * a[8 + 3*BLOCK];
    2155 
    2156       t0 -= region2[9] * a[9 + 0*BLOCK];
    2157       t1 -= region2[9] * a[9 + 1*BLOCK];
    2158       t2 -= region2[9] * a[9 + 2*BLOCK];
    2159       t3 -= region2[9] * a[9 + 3*BLOCK];
    2160 
    2161       t0 -= region2[10] * a[10 + 0*BLOCK];
    2162       t1 -= region2[10] * a[10 + 1*BLOCK];
    2163       t2 -= region2[10] * a[10 + 2*BLOCK];
    2164       t3 -= region2[10] * a[10 + 3*BLOCK];
    2165 
    2166       t0 -= region2[11] * a[11 + 0*BLOCK];
    2167       t1 -= region2[11] * a[11 + 1*BLOCK];
    2168       t2 -= region2[11] * a[11 + 2*BLOCK];
    2169       t3 -= region2[11] * a[11 + 3*BLOCK];
    2170 
    2171       t0 -= region2[12] * a[12 + 0*BLOCK];
    2172       t1 -= region2[12] * a[12 + 1*BLOCK];
    2173       t2 -= region2[12] * a[12 + 2*BLOCK];
    2174       t3 -= region2[12] * a[12 + 3*BLOCK];
    2175 
    2176       t0 -= region2[13] * a[13 + 0*BLOCK];
    2177       t1 -= region2[13] * a[13 + 1*BLOCK];
    2178       t2 -= region2[13] * a[13 + 2*BLOCK];
    2179       t3 -= region2[13] * a[13 + 3*BLOCK];
    2180 
    2181       t0 -= region2[14] * a[14 + 0*BLOCK];
    2182       t1 -= region2[14] * a[14 + 1*BLOCK];
    2183       t2 -= region2[14] * a[14 + 2*BLOCK];
    2184       t3 -= region2[14] * a[14 + 3*BLOCK];
    2185 
    2186       t0 -= region2[15] * a[15 + 0*BLOCK];
    2187       t1 -= region2[15] * a[15 + 1*BLOCK];
    2188       t2 -= region2[15] * a[15 + 2*BLOCK];
    2189       t3 -= region2[15] * a[15 + 3*BLOCK];
    2190 #endif
    2191       region[0] = t0;
    2192       region[1] = t1;
    2193       region[2] = t2;
    2194       region[3] = t3;
    2195       a+=4*BLOCK;
    2196       region+=4;
    2197     }
    2198   } else {
    2199 #endif
    2200     for (j = 0; j < BLOCK; j ++) {
    2201       CoinWorkDouble t00 = region[j];
    2202       for (k = 0; k < n; ++k) {
    2203         t00 -= region2[k] * a[k + j * BLOCK];
    2204       }
    2205       region[j] = t00;
    2206     }
    2207 #ifdef BLOCKUNROLL
    2208   }
    2209 #endif
    2210 }
  • trunk/Clp/src/ClpCholeskyDense.hpp

    r1370 r1371  
    2323  /** Factorize - filling in rowsDropped and returning number dropped.
    2424      If return code negative then out of memory */
    25   virtual int factorize(const double * diagonal, int * rowsDropped) ;
     25  virtual int factorize(const CoinWorkDouble * diagonal, int * rowsDropped) ;
    2626  /** Uses factorization to solve. */
    27   virtual void solve (double * region) ;
     27  virtual void solve (CoinWorkDouble * region) ;
    2828  /**@}*/
    2929
     
    4141  void factorizePart3(int * rowsDropped) ;
    4242  /** Forward part of solve */
    43   void solveF1(longDouble * a,int n,double * region);
    44   void solveF2(longDouble * a,int n,double * region,double * region2);
     43  void solveF1(longDouble * a,int n,CoinWorkDouble * region);
     44  void solveF2(longDouble * a,int n,CoinWorkDouble * region,CoinWorkDouble * region2);
    4545  /** Backward part of solve */
    46   void solveB1(longDouble * a,int n,double * region);
    47   void solveB2(longDouble * a,int n,double * region,double * region2);
    48   /** Uses factorization to solve. */
    49   void solveLong (CoinWorkDouble * region) ;
    50   /** Forward part of solve */
    51   void solveF1Long(longDouble * a,int n,CoinWorkDouble * region);
    52   void solveF2Long(longDouble * a,int n,CoinWorkDouble * region,CoinWorkDouble * region2);
    53   /** Backward part of solve */
    54   void solveB1Long(longDouble * a,int n,CoinWorkDouble * region);
    55   void solveB2Long(longDouble * a,int n,CoinWorkDouble * region,CoinWorkDouble * region2);
    56   /** Uses factorization to solve. */
    57   void solveLongWork (CoinWorkDouble * region) ;
    58   /** Forward part of solve */
    59   void solveF1LongWork(longDouble * a,int n,CoinWorkDouble * region);
    60   void solveF2LongWork(longDouble * a,int n,CoinWorkDouble * region,CoinWorkDouble * region2);
    61   /** Backward part of solve */
    62   void solveB1LongWork(longDouble * a,int n,CoinWorkDouble * region);
    63   void solveB2LongWork(longDouble * a,int n,CoinWorkDouble * region,CoinWorkDouble * region2);
     46  void solveB1(longDouble * a,int n,CoinWorkDouble * region);
     47  void solveB2(longDouble * a,int n,CoinWorkDouble * region,CoinWorkDouble * region2);
    6448  int bNumber(const longDouble * array,int &, int&);
    6549  /** A */
  • trunk/Clp/src/ClpFactorization.cpp

    r1370 r1371  
    11211121  coinFactorizationB_ = NULL;
    11221122  //coinFactorizationB_ = new CoinSmallFactorization();
     1123  forceB_=0;
     1124  goOslThreshold_ = -1;
    11231125  goDenseThreshold_ = -1;
    11241126  goSmallThreshold_ = -1;
     
    11401142    networkBasis_=NULL;
    11411143#endif
     1144  forceB_ = rhs.forceB_;
     1145  goOslThreshold_ = rhs.goOslThreshold_;
    11421146  goDenseThreshold_ = rhs.goDenseThreshold_;
    11431147  goSmallThreshold_ = rhs.goSmallThreshold_;
     
    11481152    else if (denseIfSmaller<=goSmallThreshold_)
    11491153      goDense=2;
     1154    else if (denseIfSmaller<=goOslThreshold_)
     1155      goDense=3;
    11501156  } else if (denseIfSmaller<0) {
    11511157    if (-denseIfSmaller<=goDenseThreshold_)
     
    11531159    else if (-denseIfSmaller<=goSmallThreshold_)
    11541160      goDense=2;
     1161    else if (-denseIfSmaller<=goOslThreshold_)
     1162      goDense=3;
    11551163  }
    11561164  if (rhs.coinFactorizationA_&&!goDense)
     
    11661174    if (goDense==1)
    11671175      coinFactorizationB_ = new CoinDenseFactorization();
     1176    else if (goDense==2)
     1177      coinFactorizationB_ = new CoinSimpFactorization();
    11681178    else
    1169       coinFactorizationB_ = new CoinSimpFactorization();
     1179      coinFactorizationB_ = new CoinOslFactorization();
    11701180    if (rhs.coinFactorizationA_) {
    11711181      coinFactorizationB_->maximumPivots(rhs.coinFactorizationA_->maximumPivots());
     
    11981208  factorization_instrument(1);
    11991209#endif
     1210  forceB_=0;
     1211  goOslThreshold_ = -1;
    12001212  goDenseThreshold_ = -1;
    12011213  goSmallThreshold_ = -1;
     
    12141226  coinFactorizationB_ = rhs.clone();
    12151227  //coinFactorizationB_ = new CoinSmallFactorization(rhs);
     1228  forceB_=0;
     1229  goOslThreshold_ = -1;
    12161230  goDenseThreshold_ = -1;
    12171231  goSmallThreshold_ = -1;
     
    12511265      networkBasis_=NULL;
    12521266#endif
    1253     delete coinFactorizationA_;
     1267    forceB_ = rhs.forceB_;
     1268    goOslThreshold_ = rhs.goOslThreshold_;
    12541269    goDenseThreshold_ = rhs.goDenseThreshold_;
    12551270    goSmallThreshold_ = rhs.goSmallThreshold_;
    1256     if (rhs.coinFactorizationA_)
    1257       coinFactorizationA_ = new CoinFactorization(*(rhs.coinFactorizationA_));
    1258     else
     1271    if (rhs.coinFactorizationA_) {
     1272      if (coinFactorizationA_)
     1273        *coinFactorizationA_ = *(rhs.coinFactorizationA_);
     1274      else
     1275        coinFactorizationA_ = new CoinFactorization(*rhs.coinFactorizationA_);
     1276    } else {
     1277      delete coinFactorizationA_;
    12591278      coinFactorizationA_=NULL;
    1260     delete coinFactorizationB_;
     1279    }
     1280   
    12611281    if (rhs.coinFactorizationB_) {
    1262       coinFactorizationB_ = rhs.coinFactorizationB_->clone();
    1263       //coinFactorizationB_ = new CoinSmallFactorization(*(rhs.coinFactorizationB_));
     1282      if (coinFactorizationB_) {
     1283        CoinDenseFactorization * denseR = dynamic_cast<CoinDenseFactorization *>(rhs.coinFactorizationB_);
     1284        CoinDenseFactorization * dense = dynamic_cast<CoinDenseFactorization *>(coinFactorizationB_);
     1285        CoinOslFactorization * oslR = dynamic_cast<CoinOslFactorization *>(rhs.coinFactorizationB_);
     1286        CoinOslFactorization * osl = dynamic_cast<CoinOslFactorization *>(coinFactorizationB_);
     1287        CoinSimpFactorization * simpR = dynamic_cast<CoinSimpFactorization *>(rhs.coinFactorizationB_);
     1288        CoinSimpFactorization * simp = dynamic_cast<CoinSimpFactorization *>(coinFactorizationB_);
     1289        if (dense&&denseR) {
     1290          *dense=*denseR;
     1291        } else if (osl&&oslR) {
     1292          *osl=*oslR;
     1293        } else if (simp&&simpR) {
     1294          *simp=*simpR;
     1295        } else {
     1296          delete coinFactorizationB_;
     1297          coinFactorizationB_ = rhs.coinFactorizationB_->clone();
     1298        }
     1299      } else {
     1300        coinFactorizationB_ = rhs.coinFactorizationB_->clone();
     1301      }
    12641302    } else {
     1303      delete coinFactorizationB_;
    12651304      coinFactorizationB_=NULL;
    12661305    }
     
    12761315ClpFactorization::goDenseOrSmall(int numberRows)
    12771316{
    1278   if (numberRows<=goDenseThreshold_) {
    1279     delete coinFactorizationA_;
    1280     delete coinFactorizationB_;
    1281     coinFactorizationA_ = NULL;
    1282     coinFactorizationB_ = new CoinDenseFactorization();
    1283     //printf("going dense\n");
    1284   } else if (numberRows<=goSmallThreshold_) {
    1285     delete coinFactorizationA_;
    1286     delete coinFactorizationB_;
    1287     coinFactorizationA_ = NULL;
    1288     coinFactorizationB_ = new CoinSimpFactorization();
    1289     //printf("going small\n");
     1317  if (!forceB_) {
     1318    if (numberRows<=goDenseThreshold_) {
     1319      delete coinFactorizationA_;
     1320      delete coinFactorizationB_;
     1321      coinFactorizationA_ = NULL;
     1322      coinFactorizationB_ = new CoinDenseFactorization();
     1323      //printf("going dense\n");
     1324    } else if (numberRows<=goSmallThreshold_) {
     1325      delete coinFactorizationA_;
     1326      delete coinFactorizationB_;
     1327      coinFactorizationA_ = NULL;
     1328      coinFactorizationB_ = new CoinSimpFactorization();
     1329      //printf("going small\n");
     1330    } else if (numberRows<=goOslThreshold_) {
     1331      delete coinFactorizationA_;
     1332      delete coinFactorizationB_;
     1333      coinFactorizationA_ = NULL;
     1334      coinFactorizationB_ = new CoinOslFactorization();
     1335      //printf("going small\n");
     1336    }
    12901337  }
    12911338  assert (!coinFactorizationA_||!coinFactorizationB_);
     1339}
     1340// If nonzero force use of 1,dense 2,small 3,osl
     1341void
     1342ClpFactorization::forceOtherFactorization(int which)
     1343{
     1344  delete coinFactorizationB_;
     1345  forceB_=0;
     1346  coinFactorizationB_ = NULL;
     1347  if (which>0&&which<4) {
     1348    delete coinFactorizationA_;
     1349    coinFactorizationA_ = NULL;
     1350    forceB_=which;
     1351    switch (which) {
     1352    case 1:
     1353      coinFactorizationB_ = new CoinDenseFactorization();
     1354      goDenseThreshold_=COIN_INT_MAX;
     1355    break;
     1356    case 2:
     1357      coinFactorizationB_ = new CoinSimpFactorization();
     1358      goSmallThreshold_=COIN_INT_MAX;
     1359    break;
     1360    case 3:
     1361      coinFactorizationB_ = new CoinOslFactorization();
     1362      goOslThreshold_=COIN_INT_MAX;
     1363    break;
     1364    }
     1365  } else if (!coinFactorizationA_) {
     1366    coinFactorizationA_ = new CoinFactorization();
     1367  }
    12921368}
    12931369int
     
    13131389    int solveMode=2;
    13141390    if (model->numberIterations())
    1315       solveMode+=100;
     1391      solveMode+=8;
    13161392    if (valuesPass)
    1317       solveMode+= 10;
    1318     coinFactorizationB_->setSolveMode(-1);
     1393      solveMode+= 4;
     1394    coinFactorizationB_->setSolveMode(solveMode);
    13191395    while (status()<-98) {
    13201396     
     
    13861462        numberBasic=numberRowBasic;
    13871463        matrix->generalExpanded(model,0,numberBasic);
     1464      } else if (numberBasic<numberRows) {
     1465        // add in some
     1466        int needed = numberRows-numberBasic;
     1467        // move up columns
     1468        for (i=numberBasic-1;i>=numberRowBasic;i--)
     1469          pivotTemp[i+needed]=pivotTemp[i];
     1470        numberRowBasic=0;
     1471        numberBasic=numberRows;
     1472        for (i=0;i<numberRows;i++) {
     1473          if (model->getRowStatus(i) == ClpSimplex::basic) {
     1474            pivotTemp[numberRowBasic++]=i;
     1475          } else if (needed) {
     1476            needed--;
     1477            model->setRowStatus(i,ClpSimplex::basic);
     1478            pivotTemp[numberRowBasic++]=i;
     1479          }
     1480        }
    13881481      }
    13891482      CoinBigIndex numberElements=numberRowBasic;
     
    13991492      // Not needed for dense
    14001493      numberElements = 3 * numberBasic + 3 * numberElements + 20000;
     1494      int numberIterations=model->numberIterations();
     1495      coinFactorizationB_->setUsefulInformation(&numberIterations,0);
    14011496      coinFactorizationB_->getAreas ( numberRows,
    14021497                                      numberRowBasic+numberColumnBasic, numberElements,
     
    14551550      coinFactorizationB_->factor (  );
    14561551      if (coinFactorizationB_->status() == -1 &&
    1457           (coinFactorizationB_->solveMode()%10)!=0) {
     1552          (coinFactorizationB_->solveMode()%3)!=0) {
    14581553        int solveMode=coinFactorizationB_->solveMode();
    1459         solveMode -= solveMode%10; // so bottom will be 0
     1554        solveMode -= solveMode%3; // so bottom will be 0
    14601555        coinFactorizationB_->setSolveMode(solveMode);
    14611556        coinFactorizationB_->setStatus(-99);
     1557      }
     1558      if (coinFactorizationB_->status() == -99)
    14621559        continue;
    1463       }
    14641560      // If we get here status is 0 or -1
    14651561      if (coinFactorizationB_->status() == 0&&numberBasic==numberRows) {
    14661562        coinFactorizationB_->postProcess(pivotTemp,pivotVariable);
    1467       } else {
     1563      } else if (solveType==0||solveType==2/*||solveType==1*/) {
    14681564        // Change pivotTemp to be correct list
    14691565        anyChanged=true;
     
    21192215  return coinFactorizationA_->status();
    21202216}
    2121 /* Replaces one Column to basis,
     2217/* Replaces one Column in basis,
    21222218   returns 0=OK, 1=Probably OK, 2=singular, 3=no room
    21232219   If checkBeforeModifying is true will do all accuracy checks
     
    21492245                                                        checkBeforeModifying);
    21502246      } else {
    2151         returnCode= coinFactorizationB_->replaceColumn(tableauColumn,
    2152                                               pivotRow,
    2153                                               pivotCheck,
    2154                                               checkBeforeModifying);
     2247        bool tab = coinFactorizationB_->wantsTableauColumn();
     2248        int numberIterations=model->numberIterations();
     2249        coinFactorizationB_->setUsefulInformation(&numberIterations,1);
     2250        returnCode=
     2251          coinFactorizationB_->replaceColumn(tab?tableauColumn:regionSparse,
     2252                                             pivotRow,
     2253                                             pivotCheck,
     2254                                             checkBeforeModifying);
    21552255#ifdef CLP_DEBUG
    21562256        // check basic
     
    23162416    factorization_instrument(5);
    23172417#endif
     2418    //#define PRINT_VECTOR
     2419#ifdef PRINT_VECTOR
     2420    printf("Update\n");
     2421    regionSparse2->print();
     2422#endif
    23182423    return returnCode;
    23192424#ifndef SLIM_CLP
     
    24272532    factorization_instrument(9);
    24282533#endif
     2534#ifdef PRINT_VECTOR
     2535    printf("UpdateTwoFT\n");
     2536    regionSparse2->print();
     2537    regionSparse3->print();
     2538#endif
    24292539    return returnCode;
    24302540#ifndef SLIM_CLP
     
    24692579#endif
    24702580    int returnCode;
     2581
    24712582    if (coinFactorizationA_) {
    24722583      coinFactorizationA_->setCollectStatistics(true);
     
    24802591#ifdef CLP_FACTORIZATION_INSTRUMENT
    24812592    factorization_instrument(6);
     2593#endif
     2594#ifdef PRINT_VECTOR
     2595    printf("UpdateTranspose\n");
     2596    regionSparse2->print();
    24822597#endif
    24832598    return returnCode;
     
    26512766  pivotTolerance(CoinMin(CoinMax(pivotTolerance(),newValue),0.999));
    26522767}
    2653 #endif
     2768// Sets factorization
     2769void
     2770ClpFactorization::setFactorization(ClpFactorization & rhs)
     2771{
     2772  ClpFactorization::operator=(rhs);
     2773}
     2774#endif
  • trunk/Clp/src/ClpFactorization.hpp

    r1370 r1371  
    1313class ClpNetworkBasis;
    1414#ifndef CLP_MULTIPLE_FACTORIZATIONS
     15#ifdef CLP_OSL
     16#define CLP_MULTIPLE_FACTORIZATIONS 4
     17#else
    1518#define CLP_MULTIPLE_FACTORIZATIONS 3
     19#endif
    1620#endif   
    1721#if CLP_MULTIPLE_FACTORIZATIONS == 1
    1822#include "CoinDenseFactorization.hpp"
    1923typedef CoinDenseFactorization CoinSmallFactorization;
     24typedef CoinOslFactorization CoinSmallFactorization;
    2025#elif CLP_MULTIPLE_FACTORIZATIONS == 2
    2126#include "CoinSimpFactorization.hpp"
    2227typedef CoinSimpFactorization CoinSmallFactorization;
     28typedef CoinOslFactorization CoinSmallFactorization;
    2329#elif CLP_MULTIPLE_FACTORIZATIONS == 3
    2430#include "CoinDenseFactorization.hpp"
    2531#include "CoinSimpFactorization.hpp"
     32#define CoinOslFactorization CoinDenseFactorization
     33#elif CLP_MULTIPLE_FACTORIZATIONS == 4
     34#include "CoinDenseFactorization.hpp"
     35#include "CoinSimpFactorization.hpp"
     36#include "CoinOslFactorization.hpp"
    2637#endif
    2738
     
    275286    }
    276287  }
     288  /// If nonzero force use of 1,dense 2,small 3,osl
     289  void forceOtherFactorization(int which);
     290  /// Get switch to osl if number rows <= this
     291  inline int goOslThreshold() const
     292  { return goOslThreshold_;}
     293  /// Set switch to osl if number rows <= this
     294  inline void setGoOslThreshold(int value)
     295  { goOslThreshold_ = value;}
    277296  /// Get switch to dense if number rows <= this
    278297  inline int goDenseThreshold() const
     
    289308  /// Go over to dense or small code if small enough
    290309  void goDenseOrSmall(int numberRows) ;
     310  /// Sets factorization
     311  void setFactorization(ClpFactorization & factorization);
    291312  /// Return 1 if dense code
    292313  inline int isDenseOrSmall() const
     
    347368  /// Pointer to CoinSmallFactorization
    348369  CoinSmallFactorization * coinFactorizationB_;
     370  /// If nonzero force use of 1,dense 2,small 3,osl
     371  int forceB_;
     372  /// Switch to osl if number rows <= this
     373  int goOslThreshold_;
    349374  /// Switch to small if number rows <= this
    350375  int goSmallThreshold_;
  • trunk/Clp/src/ClpMain.cpp

    r1370 r1371  
    520520            case PFI:
    521521              models[iModel].factorization()->setForrestTomlin(action==0);
     522              break;
     523            case FACTORIZATION:
     524              models[iModel].factorization()->forceOtherFactorization(action);
    522525              break;
    523526            case CRASH:
  • trunk/Clp/src/ClpNode.cpp

    r1370 r1371  
    9595    maximumColumns_=CoinMax(maximumColumns_,numberColumns);
    9696    maximumTotal = maximumRows_+maximumColumns_;
     97    assert (!factorization_);
    9798    factorization_ = new ClpFactorization(*model->factorization(),numberRows);
    9899    status_ = CoinCopyOfArrayPartial(model->statusArray(),maximumTotal,numberTotal);
  • trunk/Clp/src/ClpPrimalColumnSteepest.cpp

    r1370 r1371  
    15551555      for (j=0;j<number;j++) {
    15561556        int iSequence = index[j];
     1557        assert (iSequence>=0&&iSequence<model_->numberRows());
    15571558        double thisWeight = weight[iSequence];
    15581559        // row has -1
  • trunk/Clp/src/ClpSimplex.cpp

    r1370 r1371  
    22102210    saveStatus_ = NULL;
    22112211  }
    2212   delete factorization_;
    22132212  if (rhs.factorization_) {
    2214     factorization_ = new ClpFactorization(*rhs.factorization_,numberRows_);
     2213    setFactorization(*rhs.factorization_);
    22152214  } else {
     2215    delete factorization_;
    22162216    factorization_=NULL;
    22172217  }
     
    23432343    }
    23442344  }
    2345   delete rowCopy_;
    2346   rowCopy_=NULL;
    23472345  delete [] saveStatus_;
    23482346  saveStatus_=NULL;
     2347  if (type!=1) {
     2348    delete rowCopy_;
     2349    rowCopy_=NULL;
     2350  }
    23492351  if (!type) {
    23502352    // delete everything
     
    29652967      return false;
    29662968    }
    2967     if (makeRowCopy&&!oldMatrix) {
    2968       delete rowCopy_;
    2969       // may return NULL if can't give row copy
    2970       rowCopy_ = matrix_->reverseOrderedCopy();
     2969    bool rowCopyIsScaled;
     2970    if (makeRowCopy) {
     2971      if(!oldMatrix||!rowCopy_) {
     2972        delete rowCopy_;
     2973        // may return NULL if can't give row copy
     2974        rowCopy_ = matrix_->reverseOrderedCopy();
     2975        rowCopyIsScaled=false;
     2976      } else {
     2977        rowCopyIsScaled=true;
     2978      }
    29712979    }
    29722980#if 0
     
    31313139        }
    31323140      }
    3133     } else if (makeRowCopy&&scalingFlag_>0&&!oldMatrix) {
     3141    } else if (makeRowCopy&&scalingFlag_>0&&!rowCopyIsScaled) {
    31343142      matrix_->scaleRowCopy(this);
    31353143    }
     
    36883696      int iRow,iColumn;
    36893697      for (iRow=0;iRow<4;iRow++) {
    3690         rowArray_[iRow]->clear();
    3691 #ifndef NDEBUG
    36923698        int length =numberRows2+factorization_->maximumPivots();
    36933699        if (iRow==3||objective_->type()>1)
    36943700          length += numberColumns_;
    3695         assert(rowArray_[iRow]->capacity()>=length);
     3701        if(rowArray_[iRow]->capacity()>=length) {
     3702          rowArray_[iRow]->clear();
     3703        } else {
     3704          // model size or maxinv changed
     3705          rowArray_[iRow]->reserve(length);
     3706        }
     3707#ifndef NDEBUG
    36963708        rowArray_[iRow]->checkClear();
    36973709#endif
     
    36993711     
    37003712      for (iColumn=0;iColumn<2;iColumn++) {
    3701         columnArray_[iColumn]->clear();
    3702 #ifndef NDEBUG
    37033713        int length =numberColumns_;
    37043714        if (iColumn)
    37053715          length=CoinMax(numberRows2,numberColumns_);
    3706         assert(columnArray_[iColumn]->capacity()>=length);
     3716        if(columnArray_[iColumn]->capacity()>=length) {
     3717          columnArray_[iColumn]->clear();
     3718        } else {
     3719          // model size or maxinv changed
     3720          columnArray_[iColumn]->reserve(length);
     3721        }
     3722#ifndef NDEBUG
    37073723        columnArray_[iColumn]->checkClear();
    37083724#endif
     
    43104326  primalColumnPivot_ = choice.clone(true);
    43114327}
    4312 // Passes in factorization
    43134328void
    43144329ClpSimplex::setFactorization( ClpFactorization & factorization)
    43154330{
    4316   delete factorization_;
    4317   factorization_= new ClpFactorization(factorization,numberRows_);
     4331  if (factorization_)
     4332    factorization_->setFactorization(factorization);
     4333  else
     4334    factorization_ = new ClpFactorization(factorization,
     4335                                          numberRows_);
     4336}
     4337
     4338// Swaps factorization
     4339ClpFactorization *
     4340ClpSimplex::swapFactorization( ClpFactorization * factorization)
     4341{
     4342  ClpFactorization * swap =factorization_;
     4343  factorization_= factorization;
     4344  return swap;
    43184345}
    43194346// Copies in factorization to existing one
  • trunk/Clp/src/ClpSimplex.hpp

    r1370 r1371  
    364364  /// Passes in factorization
    365365  void setFactorization( ClpFactorization & factorization);
     366  // Swaps factorization
     367  ClpFactorization * swapFactorization( ClpFactorization * factorization);
    366368  /// Copies in factorization to existing one
    367369  void copyFactorization( ClpFactorization & factorization);
  • trunk/Clp/src/ClpSimplexDual.cpp

    r1370 r1371  
    358358    } else if (!ifValuesPass) {
    359359      gutsOfSolution(NULL,NULL);
     360      // double check
     361      if (numberDualInfeasibilities_||numberPrimalInfeasibilities_)
     362        problemStatus_=-1;
    360363    }
    361364    if (usePrimal) {
     
    11001103    if(numberIterations_==debugIteration) {
    11011104      printf("dodgy iteration coming up\n");
     1105    }
     1106#endif
     1107#if 0
     1108    printf("checking nz\n");
     1109    for (int i=0;i<3;i++) {
     1110      if (!rowArray_[i]->getNumElements())
     1111        rowArray_[i]->checkClear();
    11021112    }
    11031113#endif
     
    56725682    columnUpper_[iColumn] = saveBound;
    56735683    CoinMemcpyN(savePivot, numberRows_,pivotVariable_);
    5674     delete factorization_;
    5675     factorization_ = new ClpFactorization(saveFactorization,numberRows_);
    5676 
     5684    //delete factorization_;
     5685    //factorization_ = new ClpFactorization(saveFactorization,numberRows_);
     5686    setFactorization(saveFactorization);
    56775687    newUpper[i]=objectiveChange;
    56785688#ifdef CLP_DEBUG
     
    57385748    columnLower_[iColumn] = saveBound;
    57395749    CoinMemcpyN(savePivot, numberRows_,pivotVariable_);
    5740     delete factorization_;
    5741     factorization_ = new ClpFactorization(saveFactorization,numberRows_);
     5750    //delete factorization_;
     5751    //factorization_ = new ClpFactorization(saveFactorization,numberRows_);
     5752    setFactorization(saveFactorization);
    57425753
    57435754    newLower[i]=objectiveChange;
Note: See TracChangeset for help on using the changeset viewer.