11 #ifndef CoinFactorization_H
12 #define CoinFactorization_H
87 int rowIsBasic[],
int columnIsBasic[] ,
104 const int indicesRow[],
105 const int indicesColumn[],
const double elements[] ,
116 int * indicesColumn[],
133 inline int status ( )
const {
268 #ifndef COIN_FAST_CODE
372 bool checkBeforeModifying=
false,
373 double acceptablePivot=1.0e-8);
380 int internalPivotRow);
398 bool noPermute=
false)
const;
407 bool noPermuteRegion3=
false) ;
439 int indicesColumn[],
double elements[] );
444 int indicesRow[],
double elements[] );
449 int indicesColumn[],
double elements[] );
459 int replaceRow (
int whichRow,
int numberElements,
460 const int indicesColumn[],
const double elements[] );
462 void emptyRows(
int numberToEmpty,
const int which[]);
495 int possibleDuplicates = -1 );
557 inline void addLink (
int index,
int count ) {
561 int next = firstCount[count];
562 lastCount[index] = -2 - count;
565 firstCount[count] = index;
566 nextCount[index] = -1;
568 firstCount[count] = index;
569 nextCount[index] = next;
570 lastCount[next] = index;
577 int next = nextCount[index];
578 int last = lastCount[index];
580 nextCount[last] = next;
582 int count = -last - 2;
584 firstCount[count] = next;
587 lastCount[next] = last;
589 nextCount[index] = -2;
590 lastCount[index] = -2;
618 int * indexIn)
const;
621 int * indexIn)
const;
627 int & numberNonZero1,
630 int & numberNonZero2,
644 int smallestIndex)
const;
648 int smallestIndex)
const;
652 int smallestIndex)
const;
659 int smallestIndex)
const;
684 int pivotRow,
double alpha);
688 int checkPivot(
double saveFromU,
double oldPivot)
const;
691 #define COINFACTORIZATION_BITS_PER_INT 64
692 #define COINFACTORIZATION_SHIFT_PER_INT 6
693 #define COINFACTORIZATION_MASK_PER_INT 0x3f
695 #define COINFACTORIZATION_BITS_PER_INT 32
696 #define COINFACTORIZATION_SHIFT_PER_INT 5
697 #define COINFACTORIZATION_MASK_PER_INT 0x1f
699 template <
class T>
inline bool
705 unsigned int workArea2[],
724 int numberInPivotRow = numberInRow[pivotRow] - 1;
726 int numberInPivotColumn = numberInColumn[
pivotColumn] - 1;
727 CoinBigIndex endColumn = startColumn + numberInPivotColumn + 1;
732 if ( pivotColumnPosition < 0 ) {
733 for ( pivotColumnPosition = startRow; pivotColumnPosition < endRow; pivotColumnPosition++ ) {
734 int iColumn = indexColumnU[pivotColumnPosition];
735 if ( iColumn != pivotColumn ) {
736 saveColumn[put++] = iColumn;
742 for (
CoinBigIndex i = startRow ; i < pivotColumnPosition ; i++ ) {
743 saveColumn[put++] = indexColumnU[i];
746 assert (pivotColumnPosition<endRow);
747 assert (indexColumnU[pivotColumnPosition]==pivotColumn);
748 pivotColumnPosition++;
749 for ( ; pivotColumnPosition < endRow; pivotColumnPosition++ ) {
750 saveColumn[put++] = indexColumnU[pivotColumnPosition];
753 int next = nextRow[pivotRow];
754 int last = lastRow[pivotRow];
756 nextRow[last] = next;
757 lastRow[next] = last;
759 lastRow[pivotRow] = -2;
760 numberInRow[pivotRow] = 0;
767 printf(
"more memory needed in middle of invert\n");
778 if ( pivotRowPosition < 0 ) {
779 for ( pivotRowPosition = startColumn; pivotRowPosition < endColumn; pivotRowPosition++ ) {
780 int iRow = indexRowU[pivotRowPosition];
781 if ( iRow != pivotRow ) {
783 elementL[l] = elementU[pivotRowPosition];
784 markRow[iRow] =
static_cast<T
>(l - lSave);
791 while ( indexColumnU[where] != pivotColumn ) {
795 if ( where >= end ) {
799 indexColumnU[where] = indexColumnU[end - 1];
808 for ( i = startColumn; i < pivotRowPosition; i++ ) {
809 int iRow = indexRowU[i];
811 markRow[iRow] =
static_cast<T
>(l - lSave);
813 elementL[l] = elementU[i];
820 while ( indexColumnU[where] != pivotColumn ) {
824 if ( where >= end ) {
828 indexColumnU[where] = indexColumnU[end - 1];
830 assert (numberInRow[iRow]>=0);
833 assert (pivotRowPosition<endColumn);
834 assert (indexRowU[pivotRowPosition]==pivotRow);
840 for ( ; pivotRowPosition < endColumn; pivotRowPosition++ ) {
841 int iRow = indexRowU[pivotRowPosition];
843 markRow[iRow] =
static_cast<T
>(l - lSave);
845 elementL[l] = elementU[pivotRowPosition];
852 while ( indexColumnU[where] != pivotColumn ) {
856 if ( where >= end ) {
860 indexColumnU[where] = indexColumnU[end - 1];
862 assert (numberInRow[iRow]>=0);
864 markRow[pivotRow] =
static_cast<T
>(largeInteger);
868 int *indexL = &indexRowL[lSave];
874 for ( j = 0; j < numberInPivotColumn; j++ ) {
875 multipliersL[j] *= pivotMultiplier;
879 for ( iErase = 0; iErase < increment2 * numberInPivotRow;
881 workArea2[iErase] = 0;
883 CoinBigIndex added = numberInPivotRow * numberInPivotColumn;
884 unsigned int *temp2 = workArea2;
889 for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
890 int iColumn = saveColumn[jColumn];
892 CoinBigIndex endColumn = startColumn + numberInColumn[iColumn];
893 int iRow = indexRowU[startColumn];
902 int mark = markRow[iRow];
904 if ( mark == largeInteger+1 ) {
905 largest = fabs ( value );
906 positionLargest = put;
908 checkLargest =
false;
913 if ( mark != largeInteger ) {
919 temp2[word] = temp2[word] | ( 1 << bit );
922 thisPivotValue = value;
926 for ( i = startColumn + 1; i < endColumn; i++ ) {
929 int mark = markRow[iRow];
931 if ( mark == largeInteger+1 ) {
933 indexRowU[put] = iRow;
934 elementU[put] = value;
935 if ( checkLargest ) {
936 double absValue = fabs ( value );
938 if ( absValue > largest ) {
940 positionLargest = put;
944 }
else if ( mark != largeInteger ) {
950 temp2[word] = temp2[word] | ( 1 << bit );
953 thisPivotValue = value;
957 elementU[put] = elementU[startColumn];
958 indexRowU[put] = indexRowU[startColumn];
959 if ( positionLargest == startColumn ) {
960 positionLargest = put;
963 elementU[startColumn] = thisPivotValue;
964 indexRowU[startColumn] = pivotRow;
967 numberInColumn[iColumn] = put - startColumn;
969 numberInColumnPlus[iColumn]++;
970 startColumnU[iColumn]++;
972 int next = nextColumn[iColumn];
975 space = startColumnU[next] - put - numberInColumnPlus[next];
977 if ( numberInPivotColumn > space ) {
983 positionLargest = positionLargest + startColumnU[iColumn] - startColumn;
984 startColumn = startColumnU[iColumn];
985 put = startColumn + numberInColumn[iColumn];
990 for ( j = 0; j < numberInPivotColumn; j++ ) {
991 value = work[j] - thisPivotValue * multipliersL[j];
992 double absValue = fabs ( value );
994 if ( absValue > tolerance ) {
997 elementU[put] = value;
998 indexRowU[put] = indexL[j];
999 if ( absValue > largest ) {
1001 positionLargest = put;
1010 if ( temp2[word] & ( 1 << bit ) ) {
1017 while ( indexColumnU[where] != iColumn ) {
1021 if ( where >= end ) {
1025 indexColumnU[where] = indexColumnU[end - 1];
1026 numberInRow[iRow]--;
1032 temp2[word] = temp2[word] | ( 1 << bit );
1036 numberInColumn[iColumn] = put - startColumn;
1038 if ( positionLargest >= 0 ) {
1039 value = elementU[positionLargest];
1040 iRow = indexRowU[positionLargest];
1041 elementU[positionLargest] = elementU[startColumn];
1042 indexRowU[positionLargest] = indexRowU[startColumn];
1043 elementU[startColumn] = value;
1044 indexRowU[startColumn] = iRow;
1052 temp2 += increment2;
1055 unsigned int *putBase = workArea2;
1060 while ( bigLoops ) {
1064 unsigned int *putThis = putBase;
1065 int iRow = indexL[i];
1071 for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
1072 unsigned int test = *putThis;
1074 putThis += increment2;
1075 test = 1 - ( ( test >> bit ) & 1 );
1078 int next = nextRow[iRow];
1081 space = startRowU[next] - startRowU[iRow];
1082 number += numberInRow[iRow];
1083 if ( space < number ) {
1090 next = nextRow[iRow];
1091 number = numberInRow[iRow];
1093 int saveIndex = indexColumnU[startRowU[next]];
1096 for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
1097 unsigned int test = *putThis;
1099 putThis += increment2;
1100 test = 1 - ( ( test >> bit ) & 1 );
1101 indexColumnU[end] = saveColumn[jColumn];
1105 indexColumnU[startRowU[next]] = saveIndex;
1106 markRow[iRow] =
static_cast<T
>(largeInteger+1);
1107 number = end - startRowU[iRow];
1108 numberInRow[iRow] = number;
1116 for ( bit = 0; i < numberInPivotColumn; i++, bit++ ) {
1117 unsigned int *putThis = putBase;
1118 int iRow = indexL[i];
1124 for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
1125 unsigned int test = *putThis;
1127 putThis += increment2;
1128 test = 1 - ( ( test >> bit ) & 1 );
1131 int next = nextRow[iRow];
1134 space = startRowU[next] - startRowU[iRow];
1135 number += numberInRow[iRow];
1136 if ( space < number ) {
1143 next = nextRow[iRow];
1144 number = numberInRow[iRow];
1148 saveIndex = indexColumnU[startRowU[next]];
1151 for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
1152 unsigned int test = *putThis;
1154 putThis += increment2;
1155 test = 1 - ( ( test >> bit ) & 1 );
1157 indexColumnU[end] = saveColumn[jColumn];
1160 indexColumnU[startRowU[next]] = saveIndex;
1161 markRow[iRow] =
static_cast<T
>(largeInteger+1);
1162 number = end - startRowU[iRow];
1163 numberInRow[iRow] = number;
1167 markRow[pivotRow] =
static_cast<T
>(largeInteger+1);
1186 #ifndef COIN_FAST_CODE
1191 #define slackValue_ -1.0
1445 #ifdef COIN_HAS_LAPACK
1446 #define DENSE_CODE 1
1451 typedef const int cipfint;
1456 #ifdef UGLY_COIN_FACTOR_CODING
1457 #define FAC_UNSET (FAC_SET+1)
1461 CoinBigIndex startColumnThis = startColumn[iPivotColumn];
1462 CoinBigIndex endColumn = startColumnThis + numberDoColumn + 1;
1466 if ( pivotColumnPosition < 0 ) {
1467 for ( pivotColumnPosition = startRowThis; pivotColumnPosition < endRow; pivotColumnPosition++ ) {
1468 int iColumn = indexColumn[pivotColumnPosition];
1469 if ( iColumn != iPivotColumn ) {
1470 saveColumn[put++] = iColumn;
1476 for (
CoinBigIndex i = startRowThis ; i < pivotColumnPosition ; i++ ) {
1477 saveColumn[put++] = indexColumn[i];
1480 assert (pivotColumnPosition<endRow);
1481 assert (indexColumn[pivotColumnPosition]==iPivotColumn);
1482 pivotColumnPosition++;
1483 for ( ; pivotColumnPosition < endRow; pivotColumnPosition++ ) {
1484 saveColumn[put++] = indexColumn[pivotColumnPosition];
1487 int next = nextRow[iPivotRow];
1488 int last = lastRow[iPivotRow];
1490 nextRow[last] = next;
1491 lastRow[next] = last;
1492 nextRow[iPivotRow] = numberGoodU_;
1493 lastRow[iPivotRow] = -2;
1494 numberInRow[iPivotRow] = 0;
1499 if ( l + numberDoColumn > lengthAreaL_ ) {
1501 if ((messageLevel_&4)!=0)
1502 printf(
"more memory needed in middle of invert\n");
1509 startColumnL[numberGoodL_] = l;
1511 startColumnL[numberGoodL_] = l + numberDoColumn;
1512 lengthL_ += numberDoColumn;
1513 if ( pivotRowPosition < 0 ) {
1514 for ( pivotRowPosition = startColumnThis; pivotRowPosition < endColumn; pivotRowPosition++ ) {
1515 int iRow = indexRow[pivotRowPosition];
1516 if ( iRow != iPivotRow ) {
1517 indexRowL[l] = iRow;
1518 elementL[l] = element[pivotRowPosition];
1519 markRow[iRow] = l - lSave;
1526 while ( indexColumn[where] != iPivotColumn ) {
1530 if ( where >= end ) {
1534 indexColumn[where] = indexColumn[end - 1];
1535 numberInRow[iRow]--;
1543 for ( i = startColumnThis; i < pivotRowPosition; i++ ) {
1544 int iRow = indexRow[i];
1546 markRow[iRow] = l - lSave;
1547 indexRowL[l] = iRow;
1548 elementL[l] = element[i];
1555 while ( indexColumn[where] != iPivotColumn ) {
1559 if ( where >= end ) {
1563 indexColumn[where] = indexColumn[end - 1];
1564 numberInRow[iRow]--;
1565 assert (numberInRow[iRow]>=0);
1568 assert (pivotRowPosition<endColumn);
1569 assert (indexRow[pivotRowPosition]==iPivotRow);
1573 pivotRegion_.array()[numberGoodU_] = pivotMultiplier;
1575 for ( ; pivotRowPosition < endColumn; pivotRowPosition++ ) {
1576 int iRow = indexRow[pivotRowPosition];
1578 markRow[iRow] = l - lSave;
1579 indexRowL[l] = iRow;
1580 elementL[l] = element[pivotRowPosition];
1587 while ( indexColumn[where] != iPivotColumn ) {
1591 if ( where >= end ) {
1595 indexColumn[where] = indexColumn[end - 1];
1596 numberInRow[iRow]--;
1597 assert (numberInRow[iRow]>=0);
1599 markRow[iPivotRow] = FAC_SET;
1601 numberInColumn[iPivotColumn] = 0;
1603 int *indexL = &indexRowL[lSave];
1609 for ( j = 0; j < numberDoColumn; j++ ) {
1610 multipliersL[j] *= pivotMultiplier;
1614 for ( iErase = 0; iErase < increment2 * numberDoRow;
1616 workArea2[iErase] = 0;
1619 unsigned int *temp2 = workArea2;
1620 int * nextColumn = nextColumn_.array();
1624 for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) {
1625 int iColumn = saveColumn[jColumn];
1627 CoinBigIndex endColumn = startColumnThis + numberInColumn[iColumn];
1628 int iRow = indexRow[startColumnThis];
1637 int mark = markRow[iRow];
1639 if ( mark == FAC_UNSET ) {
1640 largest = fabs ( value );
1641 positionLargest = put;
1643 checkLargest =
false;
1647 checkLargest =
true;
1648 if ( mark != FAC_SET ) {
1650 workArea[mark] = value;
1654 temp2[word] = temp2[word] | ( 1 << bit );
1657 thisPivotValue = value;
1661 for ( i = startColumnThis + 1; i < endColumn; i++ ) {
1664 int mark = markRow[iRow];
1666 if ( mark == FAC_UNSET ) {
1668 indexRow[put] = iRow;
1669 element[put] = value;
1670 if ( checkLargest ) {
1671 double absValue = fabs ( value );
1673 if ( absValue > largest ) {
1675 positionLargest = put;
1679 }
else if ( mark != FAC_SET ) {
1681 workArea[mark] = value;
1685 temp2[word] = temp2[word] | ( 1 << bit );
1688 thisPivotValue = value;
1692 element[put] = element[startColumnThis];
1693 indexRow[put] = indexRow[startColumnThis];
1694 if ( positionLargest == startColumnThis ) {
1695 positionLargest = put;
1698 element[startColumnThis] = thisPivotValue;
1699 indexRow[startColumnThis] = iPivotRow;
1702 numberInColumn[iColumn] = put - startColumnThis;
1703 int * numberInColumnPlus = numberInColumnPlus_.array();
1704 numberInColumnPlus[iColumn]++;
1705 startColumn[iColumn]++;
1707 int next = nextColumn[iColumn];
1710 space = startColumn[next] - put - numberInColumnPlus[next];
1712 if ( numberDoColumn > space ) {
1714 if ( !getColumnSpace ( iColumn, numberDoColumn ) ) {
1718 positionLargest = positionLargest + startColumn[iColumn] - startColumnThis;
1719 startColumnThis = startColumn[iColumn];
1720 put = startColumnThis + numberInColumn[iColumn];
1722 double tolerance = zeroTolerance_;
1724 int *nextCount = nextCount_.array();
1725 for ( j = 0; j < numberDoColumn; j++ ) {
1726 value = workArea[j] - thisPivotValue * multipliersL[j];
1727 double absValue = fabs ( value );
1729 if ( absValue > tolerance ) {
1731 element[put] = value;
1732 indexRow[put] = indexL[j];
1733 if ( absValue > largest ) {
1735 positionLargest = put;
1744 if ( temp2[word] & ( 1 << bit ) ) {
1751 while ( indexColumn[where] != iColumn ) {
1755 if ( where >= end ) {
1759 indexColumn[where] = indexColumn[end - 1];
1760 numberInRow[iRow]--;
1766 temp2[word] = temp2[word] | ( 1 << bit );
1770 numberInColumn[iColumn] = put - startColumnThis;
1772 if ( positionLargest >= 0 ) {
1773 value = element[positionLargest];
1774 iRow = indexRow[positionLargest];
1775 element[positionLargest] = element[startColumnThis];
1776 indexRow[positionLargest] = indexRow[startColumnThis];
1777 element[startColumnThis] = value;
1778 indexRow[startColumnThis] = iRow;
1781 if ( nextCount[iColumn + numberRows_] != -2 ) {
1783 deleteLink ( iColumn + numberRows_ );
1784 addLink ( iColumn + numberRows_, numberInColumn[iColumn] );
1786 temp2 += increment2;
1789 unsigned int *putBase = workArea2;
1794 while ( bigLoops ) {
1798 unsigned int *putThis = putBase;
1799 int iRow = indexL[i];
1805 for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) {
1806 unsigned int test = *putThis;
1808 putThis += increment2;
1809 test = 1 - ( ( test >> bit ) & 1 );
1812 int next = nextRow[iRow];
1815 space = startRow[next] - startRow[iRow];
1816 number += numberInRow[iRow];
1817 if ( space < number ) {
1818 if ( !getRowSpace ( iRow, number ) ) {
1824 next = nextRow[iRow];
1825 number = numberInRow[iRow];
1827 int saveIndex = indexColumn[startRow[next]];
1830 for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) {
1831 unsigned int test = *putThis;
1833 putThis += increment2;
1834 test = 1 - ( ( test >> bit ) & 1 );
1835 indexColumn[end] = saveColumn[jColumn];
1839 indexColumn[startRow[next]] = saveIndex;
1840 markRow[iRow] = FAC_UNSET;
1841 number = end - startRow[iRow];
1842 numberInRow[iRow] = number;
1843 deleteLink ( iRow );
1844 addLink ( iRow, number );
1850 for ( bit = 0; i < numberDoColumn; i++, bit++ ) {
1851 unsigned int *putThis = putBase;
1852 int iRow = indexL[i];
1858 for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) {
1859 unsigned int test = *putThis;
1861 putThis += increment2;
1862 test = 1 - ( ( test >> bit ) & 1 );
1865 int next = nextRow[iRow];
1868 space = startRow[next] - startRow[iRow];
1869 number += numberInRow[iRow];
1870 if ( space < number ) {
1871 if ( !getRowSpace ( iRow, number ) ) {
1877 next = nextRow[iRow];
1878 number = numberInRow[iRow];
1882 saveIndex = indexColumn[startRow[next]];
1885 for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) {
1886 unsigned int test = *putThis;
1888 putThis += increment2;
1889 test = 1 - ( ( test >> bit ) & 1 );
1891 indexColumn[end] = saveColumn[jColumn];
1894 indexColumn[startRow[next]] = saveIndex;
1895 markRow[iRow] = FAC_UNSET;
1896 number = end - startRow[iRow];
1897 numberInRow[iRow] = number;
1898 deleteLink ( iRow );
1899 addLink ( iRow, number );
1901 markRow[iPivotRow] = FAC_UNSET;
1903 deleteLink ( iPivotRow );
1904 deleteLink ( iPivotColumn + numberRows_ );
1905 totalElements_ += added;