12 #ifndef CoinFactorization_H
13 #define CoinFactorization_H
90 int rowIsBasic[],
int columnIsBasic[] ,
107 const int indicesRow[],
108 const int indicesColumn[],
const double elements[] ,
119 int * indicesColumn[],
136 inline int status ( )
const {
271 #ifndef COIN_FAST_CODE
375 bool checkBeforeModifying=
false,
376 double acceptablePivot=1.0e-8);
383 int internalPivotRow);
401 bool noPermute=
false)
const;
410 bool noPermuteRegion3=
false) ;
442 int indicesColumn[],
double elements[] );
447 int indicesRow[],
double elements[] );
452 int indicesColumn[],
double elements[] );
462 int replaceRow (
int whichRow,
int numberElements,
463 const int indicesColumn[],
const double elements[] );
465 void emptyRows(
int numberToEmpty,
const int which[]);
498 int possibleDuplicates = -1 );
560 inline void addLink (
int index,
int count ) {
564 int next = firstCount[count];
565 lastCount[index] = -2 - count;
568 firstCount[count] = index;
569 nextCount[index] = -1;
571 firstCount[count] = index;
572 nextCount[index] = next;
573 lastCount[next] = index;
580 int next = nextCount[index];
581 int last = lastCount[index];
583 nextCount[last] = next;
585 int count = -last - 2;
587 firstCount[count] = next;
590 lastCount[next] = last;
592 nextCount[index] = -2;
593 lastCount[index] = -2;
621 int * indexIn)
const;
624 int * indexIn)
const;
630 int & numberNonZero1,
633 int & numberNonZero2,
647 int smallestIndex)
const;
651 int smallestIndex)
const;
655 int smallestIndex)
const;
662 int smallestIndex)
const;
687 int pivotRow,
double alpha);
691 int checkPivot(
double saveFromU,
double oldPivot)
const;
694 #define COINFACTORIZATION_BITS_PER_INT 64
695 #define COINFACTORIZATION_SHIFT_PER_INT 6
696 #define COINFACTORIZATION_MASK_PER_INT 0x3f
698 #define COINFACTORIZATION_BITS_PER_INT 32
699 #define COINFACTORIZATION_SHIFT_PER_INT 5
700 #define COINFACTORIZATION_MASK_PER_INT 0x1f
702 template <
class T>
inline bool
708 unsigned int workArea2[],
727 int numberInPivotRow = numberInRow[pivotRow] - 1;
729 int numberInPivotColumn = numberInColumn[
pivotColumn] - 1;
730 CoinBigIndex endColumn = startColumn + numberInPivotColumn + 1;
735 if ( pivotColumnPosition < 0 ) {
736 for ( pivotColumnPosition = startRow; pivotColumnPosition < endRow; pivotColumnPosition++ ) {
737 int iColumn = indexColumnU[pivotColumnPosition];
738 if ( iColumn != pivotColumn ) {
739 saveColumn[put++] = iColumn;
745 for (
CoinBigIndex i = startRow ; i < pivotColumnPosition ; i++ ) {
746 saveColumn[put++] = indexColumnU[i];
749 assert (pivotColumnPosition<endRow);
750 assert (indexColumnU[pivotColumnPosition]==pivotColumn);
751 pivotColumnPosition++;
752 for ( ; pivotColumnPosition < endRow; pivotColumnPosition++ ) {
753 saveColumn[put++] = indexColumnU[pivotColumnPosition];
756 int next = nextRow[pivotRow];
757 int last = lastRow[pivotRow];
759 nextRow[last] = next;
760 lastRow[next] = last;
762 lastRow[pivotRow] = -2;
763 numberInRow[pivotRow] = 0;
770 printf(
"more memory needed in middle of invert\n");
781 if ( pivotRowPosition < 0 ) {
782 for ( pivotRowPosition = startColumn; pivotRowPosition < endColumn; pivotRowPosition++ ) {
783 int iRow = indexRowU[pivotRowPosition];
784 if ( iRow != pivotRow ) {
786 elementL[l] = elementU[pivotRowPosition];
787 markRow[iRow] =
static_cast<T
>(l - lSave);
794 while ( indexColumnU[where] != pivotColumn ) {
798 if ( where >= end ) {
802 indexColumnU[where] = indexColumnU[end - 1];
811 for ( i = startColumn; i < pivotRowPosition; i++ ) {
812 int iRow = indexRowU[i];
814 markRow[iRow] =
static_cast<T
>(l - lSave);
816 elementL[l] = elementU[i];
823 while ( indexColumnU[where] != pivotColumn ) {
827 if ( where >= end ) {
831 indexColumnU[where] = indexColumnU[end - 1];
833 assert (numberInRow[iRow]>=0);
836 assert (pivotRowPosition<endColumn);
837 assert (indexRowU[pivotRowPosition]==pivotRow);
843 for ( ; pivotRowPosition < endColumn; pivotRowPosition++ ) {
844 int iRow = indexRowU[pivotRowPosition];
846 markRow[iRow] =
static_cast<T
>(l - lSave);
848 elementL[l] = elementU[pivotRowPosition];
855 while ( indexColumnU[where] != pivotColumn ) {
859 if ( where >= end ) {
863 indexColumnU[where] = indexColumnU[end - 1];
865 assert (numberInRow[iRow]>=0);
867 markRow[pivotRow] =
static_cast<T
>(largeInteger);
871 int *indexL = &indexRowL[lSave];
877 for ( j = 0; j < numberInPivotColumn; j++ ) {
878 multipliersL[j] *= pivotMultiplier;
882 for ( iErase = 0; iErase < increment2 * numberInPivotRow;
884 workArea2[iErase] = 0;
886 CoinBigIndex added = numberInPivotRow * numberInPivotColumn;
887 unsigned int *temp2 = workArea2;
892 for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
893 int iColumn = saveColumn[jColumn];
895 CoinBigIndex endColumn = startColumn + numberInColumn[iColumn];
896 int iRow = indexRowU[startColumn];
905 int mark = markRow[iRow];
907 if ( mark == largeInteger+1 ) {
908 largest = fabs ( value );
909 positionLargest = put;
911 checkLargest =
false;
916 if ( mark != largeInteger ) {
922 temp2[word] = temp2[word] | ( 1 << bit );
925 thisPivotValue = value;
929 for ( i = startColumn + 1; i < endColumn; i++ ) {
932 int mark = markRow[iRow];
934 if ( mark == largeInteger+1 ) {
936 indexRowU[put] = iRow;
937 elementU[put] = value;
938 if ( checkLargest ) {
939 double absValue = fabs ( value );
941 if ( absValue > largest ) {
943 positionLargest = put;
947 }
else if ( mark != largeInteger ) {
953 temp2[word] = temp2[word] | ( 1 << bit );
956 thisPivotValue = value;
960 elementU[put] = elementU[startColumn];
961 indexRowU[put] = indexRowU[startColumn];
962 if ( positionLargest == startColumn ) {
963 positionLargest = put;
966 elementU[startColumn] = thisPivotValue;
967 indexRowU[startColumn] = pivotRow;
970 numberInColumn[iColumn] = put - startColumn;
972 numberInColumnPlus[iColumn]++;
973 startColumnU[iColumn]++;
975 int next = nextColumn[iColumn];
978 space = startColumnU[next] - put - numberInColumnPlus[next];
980 if ( numberInPivotColumn > space ) {
986 if (positionLargest >= 0)
987 positionLargest = positionLargest + startColumnU[iColumn] - startColumn;
988 startColumn = startColumnU[iColumn];
989 put = startColumn + numberInColumn[iColumn];
994 for ( j = 0; j < numberInPivotColumn; j++ ) {
995 value = work[j] - thisPivotValue * multipliersL[j];
996 double absValue = fabs ( value );
998 if ( absValue > tolerance ) {
1001 elementU[put] = value;
1002 indexRowU[put] = indexL[j];
1003 if ( absValue > largest ) {
1005 positionLargest = put;
1014 if ( temp2[word] & ( 1 << bit ) ) {
1021 while ( indexColumnU[where] != iColumn ) {
1025 if ( where >= end ) {
1029 indexColumnU[where] = indexColumnU[end - 1];
1030 numberInRow[iRow]--;
1036 temp2[word] = temp2[word] | ( 1 << bit );
1040 numberInColumn[iColumn] = put - startColumn;
1042 if ( positionLargest >= 0 ) {
1043 value = elementU[positionLargest];
1044 iRow = indexRowU[positionLargest];
1045 elementU[positionLargest] = elementU[startColumn];
1046 indexRowU[positionLargest] = indexRowU[startColumn];
1047 elementU[startColumn] = value;
1048 indexRowU[startColumn] = iRow;
1056 temp2 += increment2;
1059 unsigned int *putBase = workArea2;
1064 while ( bigLoops ) {
1068 unsigned int *putThis = putBase;
1069 int iRow = indexL[i];
1075 for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
1076 unsigned int test = *putThis;
1078 putThis += increment2;
1079 test = 1 - ( ( test >> bit ) & 1 );
1082 int next = nextRow[iRow];
1085 space = startRowU[next] - startRowU[iRow];
1086 number += numberInRow[iRow];
1087 if ( space < number ) {
1094 next = nextRow[iRow];
1095 number = numberInRow[iRow];
1097 int saveIndex = indexColumnU[startRowU[next]];
1100 for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
1101 unsigned int test = *putThis;
1103 putThis += increment2;
1104 test = 1 - ( ( test >> bit ) & 1 );
1105 indexColumnU[end] = saveColumn[jColumn];
1109 indexColumnU[startRowU[next]] = saveIndex;
1110 markRow[iRow] =
static_cast<T
>(largeInteger+1);
1111 number = end - startRowU[iRow];
1112 numberInRow[iRow] = number;
1120 for ( bit = 0; i < numberInPivotColumn; i++, bit++ ) {
1121 unsigned int *putThis = putBase;
1122 int iRow = indexL[i];
1128 for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
1129 unsigned int test = *putThis;
1131 putThis += increment2;
1132 test = 1 - ( ( test >> bit ) & 1 );
1135 int next = nextRow[iRow];
1138 space = startRowU[next] - startRowU[iRow];
1139 number += numberInRow[iRow];
1140 if ( space < number ) {
1147 next = nextRow[iRow];
1148 number = numberInRow[iRow];
1152 saveIndex = indexColumnU[startRowU[next]];
1155 for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
1156 unsigned int test = *putThis;
1158 putThis += increment2;
1159 test = 1 - ( ( test >> bit ) & 1 );
1161 indexColumnU[end] = saveColumn[jColumn];
1164 indexColumnU[startRowU[next]] = saveIndex;
1165 markRow[iRow] =
static_cast<T
>(largeInteger+1);
1166 number = end - startRowU[iRow];
1167 numberInRow[iRow] = number;
1171 markRow[pivotRow] =
static_cast<T
>(largeInteger+1);
1190 #ifndef COIN_FAST_CODE
1195 #define slackValue_ -1.0
1449 #ifdef COIN_HAS_LAPACK
1450 #define DENSE_CODE 1
1455 typedef const int cipfint;
1460 #ifdef UGLY_COIN_FACTOR_CODING
1461 #define FAC_UNSET (FAC_SET+1)
1465 CoinBigIndex startColumnThis = startColumn[iPivotColumn];
1466 CoinBigIndex endColumn = startColumnThis + numberDoColumn + 1;
1470 if ( pivotColumnPosition < 0 ) {
1471 for ( pivotColumnPosition = startRowThis; pivotColumnPosition < endRow; pivotColumnPosition++ ) {
1472 int iColumn = indexColumn[pivotColumnPosition];
1473 if ( iColumn != iPivotColumn ) {
1474 saveColumn[put++] = iColumn;
1480 for (
CoinBigIndex i = startRowThis ; i < pivotColumnPosition ; i++ ) {
1481 saveColumn[put++] = indexColumn[i];
1484 assert (pivotColumnPosition<endRow);
1485 assert (indexColumn[pivotColumnPosition]==iPivotColumn);
1486 pivotColumnPosition++;
1487 for ( ; pivotColumnPosition < endRow; pivotColumnPosition++ ) {
1488 saveColumn[put++] = indexColumn[pivotColumnPosition];
1491 int next = nextRow[iPivotRow];
1492 int last = lastRow[iPivotRow];
1494 nextRow[last] = next;
1495 lastRow[next] = last;
1496 nextRow[iPivotRow] = numberGoodU_;
1497 lastRow[iPivotRow] = -2;
1498 numberInRow[iPivotRow] = 0;
1503 if ( l + numberDoColumn > lengthAreaL_ ) {
1505 if ((messageLevel_&4)!=0)
1506 printf(
"more memory needed in middle of invert\n");
1513 startColumnL[numberGoodL_] = l;
1515 startColumnL[numberGoodL_] = l + numberDoColumn;
1516 lengthL_ += numberDoColumn;
1517 if ( pivotRowPosition < 0 ) {
1518 for ( pivotRowPosition = startColumnThis; pivotRowPosition < endColumn; pivotRowPosition++ ) {
1519 int iRow = indexRow[pivotRowPosition];
1520 if ( iRow != iPivotRow ) {
1521 indexRowL[l] = iRow;
1522 elementL[l] = element[pivotRowPosition];
1523 markRow[iRow] = l - lSave;
1530 while ( indexColumn[where] != iPivotColumn ) {
1534 if ( where >= end ) {
1538 indexColumn[where] = indexColumn[end - 1];
1539 numberInRow[iRow]--;
1547 for ( i = startColumnThis; i < pivotRowPosition; i++ ) {
1548 int iRow = indexRow[i];
1550 markRow[iRow] = l - lSave;
1551 indexRowL[l] = iRow;
1552 elementL[l] = element[i];
1559 while ( indexColumn[where] != iPivotColumn ) {
1563 if ( where >= end ) {
1567 indexColumn[where] = indexColumn[end - 1];
1568 numberInRow[iRow]--;
1569 assert (numberInRow[iRow]>=0);
1572 assert (pivotRowPosition<endColumn);
1573 assert (indexRow[pivotRowPosition]==iPivotRow);
1577 pivotRegion_.array()[numberGoodU_] = pivotMultiplier;
1579 for ( ; pivotRowPosition < endColumn; pivotRowPosition++ ) {
1580 int iRow = indexRow[pivotRowPosition];
1582 markRow[iRow] = l - lSave;
1583 indexRowL[l] = iRow;
1584 elementL[l] = element[pivotRowPosition];
1591 while ( indexColumn[where] != iPivotColumn ) {
1595 if ( where >= end ) {
1599 indexColumn[where] = indexColumn[end - 1];
1600 numberInRow[iRow]--;
1601 assert (numberInRow[iRow]>=0);
1603 markRow[iPivotRow] = FAC_SET;
1605 numberInColumn[iPivotColumn] = 0;
1607 int *indexL = &indexRowL[lSave];
1613 for ( j = 0; j < numberDoColumn; j++ ) {
1614 multipliersL[j] *= pivotMultiplier;
1618 for ( iErase = 0; iErase < increment2 * numberDoRow;
1620 workArea2[iErase] = 0;
1623 unsigned int *temp2 = workArea2;
1624 int * nextColumn = nextColumn_.array();
1628 for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) {
1629 int iColumn = saveColumn[jColumn];
1631 CoinBigIndex endColumn = startColumnThis + numberInColumn[iColumn];
1632 int iRow = indexRow[startColumnThis];
1641 int mark = markRow[iRow];
1643 if ( mark == FAC_UNSET ) {
1644 largest = fabs ( value );
1645 positionLargest = put;
1647 checkLargest =
false;
1651 checkLargest =
true;
1652 if ( mark != FAC_SET ) {
1654 workArea[mark] = value;
1658 temp2[word] = temp2[word] | ( 1 << bit );
1661 thisPivotValue = value;
1665 for ( i = startColumnThis + 1; i < endColumn; i++ ) {
1668 int mark = markRow[iRow];
1670 if ( mark == FAC_UNSET ) {
1672 indexRow[put] = iRow;
1673 element[put] = value;
1674 if ( checkLargest ) {
1675 double absValue = fabs ( value );
1677 if ( absValue > largest ) {
1679 positionLargest = put;
1683 }
else if ( mark != FAC_SET ) {
1685 workArea[mark] = value;
1689 temp2[word] = temp2[word] | ( 1 << bit );
1692 thisPivotValue = value;
1696 element[put] = element[startColumnThis];
1697 indexRow[put] = indexRow[startColumnThis];
1698 if ( positionLargest == startColumnThis ) {
1699 positionLargest = put;
1702 element[startColumnThis] = thisPivotValue;
1703 indexRow[startColumnThis] = iPivotRow;
1706 numberInColumn[iColumn] = put - startColumnThis;
1707 int * numberInColumnPlus = numberInColumnPlus_.array();
1708 numberInColumnPlus[iColumn]++;
1709 startColumn[iColumn]++;
1711 int next = nextColumn[iColumn];
1714 space = startColumn[next] - put - numberInColumnPlus[next];
1716 if ( numberDoColumn > space ) {
1718 if ( !getColumnSpace ( iColumn, numberDoColumn ) ) {
1722 positionLargest = positionLargest + startColumn[iColumn] - startColumnThis;
1723 startColumnThis = startColumn[iColumn];
1724 put = startColumnThis + numberInColumn[iColumn];
1726 double tolerance = zeroTolerance_;
1728 int *nextCount = nextCount_.array();
1729 for ( j = 0; j < numberDoColumn; j++ ) {
1730 value = workArea[j] - thisPivotValue * multipliersL[j];
1731 double absValue = fabs ( value );
1733 if ( absValue > tolerance ) {
1735 element[put] = value;
1736 indexRow[put] = indexL[j];
1737 if ( absValue > largest ) {
1739 positionLargest = put;
1748 if ( temp2[word] & ( 1 << bit ) ) {
1755 while ( indexColumn[where] != iColumn ) {
1759 if ( where >= end ) {
1763 indexColumn[where] = indexColumn[end - 1];
1764 numberInRow[iRow]--;
1770 temp2[word] = temp2[word] | ( 1 << bit );
1774 numberInColumn[iColumn] = put - startColumnThis;
1776 if ( positionLargest >= 0 ) {
1777 value = element[positionLargest];
1778 iRow = indexRow[positionLargest];
1779 element[positionLargest] = element[startColumnThis];
1780 indexRow[positionLargest] = indexRow[startColumnThis];
1781 element[startColumnThis] = value;
1782 indexRow[startColumnThis] = iRow;
1785 if ( nextCount[iColumn + numberRows_] != -2 ) {
1787 deleteLink ( iColumn + numberRows_ );
1788 addLink ( iColumn + numberRows_, numberInColumn[iColumn] );
1790 temp2 += increment2;
1793 unsigned int *putBase = workArea2;
1798 while ( bigLoops ) {
1802 unsigned int *putThis = putBase;
1803 int iRow = indexL[i];
1809 for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) {
1810 unsigned int test = *putThis;
1812 putThis += increment2;
1813 test = 1 - ( ( test >> bit ) & 1 );
1816 int next = nextRow[iRow];
1819 space = startRow[next] - startRow[iRow];
1820 number += numberInRow[iRow];
1821 if ( space < number ) {
1822 if ( !getRowSpace ( iRow, number ) ) {
1828 next = nextRow[iRow];
1829 number = numberInRow[iRow];
1831 int saveIndex = indexColumn[startRow[next]];
1834 for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) {
1835 unsigned int test = *putThis;
1837 putThis += increment2;
1838 test = 1 - ( ( test >> bit ) & 1 );
1839 indexColumn[end] = saveColumn[jColumn];
1843 indexColumn[startRow[next]] = saveIndex;
1844 markRow[iRow] = FAC_UNSET;
1845 number = end - startRow[iRow];
1846 numberInRow[iRow] = number;
1847 deleteLink ( iRow );
1848 addLink ( iRow, number );
1854 for ( bit = 0; i < numberDoColumn; i++, bit++ ) {
1855 unsigned int *putThis = putBase;
1856 int iRow = indexL[i];
1862 for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) {
1863 unsigned int test = *putThis;
1865 putThis += increment2;
1866 test = 1 - ( ( test >> bit ) & 1 );
1869 int next = nextRow[iRow];
1872 space = startRow[next] - startRow[iRow];
1873 number += numberInRow[iRow];
1874 if ( space < number ) {
1875 if ( !getRowSpace ( iRow, number ) ) {
1881 next = nextRow[iRow];
1882 number = numberInRow[iRow];
1886 saveIndex = indexColumn[startRow[next]];
1889 for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) {
1890 unsigned int test = *putThis;
1892 putThis += increment2;
1893 test = 1 - ( ( test >> bit ) & 1 );
1895 indexColumn[end] = saveColumn[jColumn];
1898 indexColumn[startRow[next]] = saveIndex;
1899 markRow[iRow] = FAC_UNSET;
1900 number = end - startRow[iRow];
1901 numberInRow[iRow] = number;
1902 deleteLink ( iRow );
1903 addLink ( iRow, number );
1905 markRow[iPivotRow] = FAC_UNSET;
1907 deleteLink ( iPivotRow );
1908 deleteLink ( iPivotColumn + numberRows_ );
1909 totalElements_ += added;