/* $Id: ClpCholeskyDense.cpp 1910 2013-01-27 02:00:13Z stefan $ */ /* Copyright (C) 2003, International Business Machines Corporation and others. All Rights Reserved. This code is licensed under the terms of the Eclipse Public License (EPL). */ #include "CoinPragma.hpp" #include "CoinHelperFunctions.hpp" #include "ClpConfig.h" #include "ClpHelperFunctions.hpp" #include "ClpInterior.hpp" #include "ClpCholeskyDense.hpp" #include "ClpMessage.hpp" #include "ClpQuadraticObjective.hpp" #if CLP_HAS_ABC #include "CoinAbcCommon.hpp" #endif #if CLP_HAS_ABC<3 #undef cilk_for #undef cilk_spawn #undef cilk_sync #define cilk_for for #define cilk_spawn #define cilk_sync #endif /*#############################################################################*/ /* Constructors / Destructor / Assignment*/ /*#############################################################################*/ /*-------------------------------------------------------------------*/ /* Default Constructor */ /*-------------------------------------------------------------------*/ ClpCholeskyDense::ClpCholeskyDense () : ClpCholeskyBase(), borrowSpace_(false) { type_ = 11;; } /*-------------------------------------------------------------------*/ /* Copy constructor */ /*-------------------------------------------------------------------*/ ClpCholeskyDense::ClpCholeskyDense (const ClpCholeskyDense & rhs) : ClpCholeskyBase(rhs), borrowSpace_(rhs.borrowSpace_) { assert(!rhs.borrowSpace_ || !rhs.sizeFactor_); /* can't do if borrowing space*/ } /*-------------------------------------------------------------------*/ /* Destructor */ /*-------------------------------------------------------------------*/ ClpCholeskyDense::~ClpCholeskyDense () { if (borrowSpace_) { /* set NULL*/ sparseFactor_ = NULL; workDouble_ = NULL; diagonal_ = NULL; } } /*----------------------------------------------------------------*/ /* Assignment operator */ /*-------------------------------------------------------------------*/ ClpCholeskyDense & ClpCholeskyDense::operator=(const ClpCholeskyDense& rhs) { if (this != &rhs) { assert(!rhs.borrowSpace_ || !rhs.sizeFactor_); /* can't do if borrowing space*/ ClpCholeskyBase::operator=(rhs); borrowSpace_ = rhs.borrowSpace_; } return *this; } /*-------------------------------------------------------------------*/ /* Clone*/ /*-------------------------------------------------------------------*/ ClpCholeskyBase * ClpCholeskyDense::clone() const { return new ClpCholeskyDense(*this); } /* If not power of 2 then need to redo a bit*/ #define BLOCK 16 #define BLOCKSHIFT 4 /* Block unroll if power of 2 and at least 8*/ #define BLOCKUNROLL #define BLOCKSQ ( BLOCK*BLOCK ) #define BLOCKSQSHIFT ( BLOCKSHIFT+BLOCKSHIFT ) #define number_blocks(x) (((x)+BLOCK-1)>>BLOCKSHIFT) #define number_rows(x) ((x)<> BLOCKSHIFT; /* allow one stripe extra*/ numberBlocks = numberBlocks + ((numberBlocks * (numberBlocks + 1)) / 2); sizeFactor_ = numberBlocks * BLOCKSQ; /*#define CHOL_COMPARE*/ #ifdef CHOL_COMPARE sizeFactor_ += 95000; #endif if (!factor) { sparseFactor_ = new longDouble [sizeFactor_]; rowsDropped_ = new char [numberRows_]; memset(rowsDropped_, 0, numberRows_); workDouble_ = new longDouble[numberRows_]; diagonal_ = new longDouble[numberRows_]; } else { borrowSpace_ = true; int numberFull = factor->numberRows(); sparseFactor_ = factor->sparseFactor() + (factor->size() - sizeFactor_); workDouble_ = factor->workDouble() + (numberFull - numberRows_); diagonal_ = factor->diagonal() + (numberFull - numberRows_); } numberRowsDropped_ = 0; return 0; } /* Returns space needed */ CoinBigIndex ClpCholeskyDense::space( int numberRows) const { int numberBlocks = (numberRows + BLOCK - 1) >> BLOCKSHIFT; /* allow one stripe extra*/ numberBlocks = numberBlocks + ((numberBlocks * (numberBlocks + 1)) / 2); CoinBigIndex sizeFactor = numberBlocks * BLOCKSQ; #ifdef CHOL_COMPARE sizeFactor += 95000; #endif return sizeFactor; } /* Orders rows and saves pointer to matrix.and model */ int ClpCholeskyDense::order(ClpInterior * model) { model_ = model; int numberRows; int numberRowsModel = model_->numberRows(); int numberColumns = model_->numberColumns(); if (!doKKT_) { numberRows = numberRowsModel; } else { numberRows = 2 * numberRowsModel + numberColumns; } reserveSpace(NULL, numberRows); rowCopy_ = model->clpMatrix()->reverseOrderedCopy(); return 0; } /* Does Symbolic factorization given permutation. This is called immediately after order. If user provides this then user must provide factorize and solve. Otherwise the default factorization is used returns non-zero if not enough memory */ int ClpCholeskyDense::symbolic() { return 0; } /* Factorize - filling in rowsDropped and returning number dropped */ int ClpCholeskyDense::factorize(const CoinWorkDouble * diagonal, int * rowsDropped) { assert (!borrowSpace_); const CoinBigIndex * columnStart = model_->clpMatrix()->getVectorStarts(); const int * columnLength = model_->clpMatrix()->getVectorLengths(); const int * row = model_->clpMatrix()->getIndices(); const double * element = model_->clpMatrix()->getElements(); const CoinBigIndex * rowStart = rowCopy_->getVectorStarts(); const int * rowLength = rowCopy_->getVectorLengths(); const int * column = rowCopy_->getIndices(); const double * elementByRow = rowCopy_->getElements(); int numberColumns = model_->clpMatrix()->getNumCols(); CoinZeroN(sparseFactor_, sizeFactor_); /*perturbation*/ CoinWorkDouble perturbation = model_->diagonalPerturbation() * model_->diagonalNorm(); perturbation = perturbation * perturbation; if (perturbation > 1.0) { #ifdef COIN_DEVELOP /*if (model_->model()->logLevel()&4) */ std::cout << "large perturbation " << perturbation << std::endl; #endif perturbation = CoinSqrt(perturbation);; perturbation = 1.0; } int iRow; int newDropped = 0; CoinWorkDouble largest = 1.0; CoinWorkDouble smallest = COIN_DBL_MAX; CoinWorkDouble delta2 = model_->delta(); /* add delta*delta to diagonal*/ delta2 *= delta2; if (!doKKT_) { longDouble * work = sparseFactor_; work--; /* skip diagonal*/ int addOffset = numberRows_ - 1; const CoinWorkDouble * diagonalSlack = diagonal + numberColumns; /* largest in initial matrix*/ CoinWorkDouble largest2 = 1.0e-20; for (iRow = 0; iRow < numberRows_; iRow++) { if (!rowsDropped_[iRow]) { CoinBigIndex startRow = rowStart[iRow]; CoinBigIndex endRow = rowStart[iRow] + rowLength[iRow]; CoinWorkDouble diagonalValue = diagonalSlack[iRow] + delta2; for (CoinBigIndex k = startRow; k < endRow; k++) { int iColumn = column[k]; CoinBigIndex start = columnStart[iColumn]; CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn]; CoinWorkDouble multiplier = diagonal[iColumn] * elementByRow[k]; for (CoinBigIndex j = start; j < end; j++) { int jRow = row[j]; if (!rowsDropped_[jRow]) { if (jRow > iRow) { CoinWorkDouble value = element[j] * multiplier; work[jRow] += value; } else if (jRow == iRow) { CoinWorkDouble value = element[j] * multiplier; diagonalValue += value; } } } } for (int j = iRow + 1; j < numberRows_; j++) largest2 = CoinMax(largest2, CoinAbs(work[j])); diagonal_[iRow] = diagonalValue; largest2 = CoinMax(largest2, CoinAbs(diagonalValue)); } else { /* dropped*/ diagonal_[iRow] = 1.0; } addOffset--; work += addOffset; } /*check sizes*/ largest2 *= 1.0e-20; largest = CoinMin(largest2, CHOL_SMALL_VALUE); int numberDroppedBefore = 0; for (iRow = 0; iRow < numberRows_; iRow++) { int dropped = rowsDropped_[iRow]; /* Move to int array*/ rowsDropped[iRow] = dropped; if (!dropped) { CoinWorkDouble diagonal = diagonal_[iRow]; if (diagonal > largest2) { diagonal_[iRow] = diagonal + perturbation; } else { diagonal_[iRow] = diagonal + perturbation; rowsDropped[iRow] = 2; numberDroppedBefore++; } } } doubleParameters_[10] = CoinMax(1.0e-20, largest); integerParameters_[20] = 0; doubleParameters_[3] = 0.0; doubleParameters_[4] = COIN_DBL_MAX; integerParameters_[34] = 0; /* say all must be positive*/ #ifdef CHOL_COMPARE if (numberRows_ < 200) factorizePart3(rowsDropped); #endif factorizePart2(rowsDropped); newDropped = integerParameters_[20] + numberDroppedBefore; largest = doubleParameters_[3]; smallest = doubleParameters_[4]; if (model_->messageHandler()->logLevel() > 1) std::cout << "Cholesky - largest " << largest << " smallest " << smallest << std::endl; choleskyCondition_ = largest / smallest; /*drop fresh makes some formADAT easier*/ if (newDropped || numberRowsDropped_) { newDropped = 0; for (int i = 0; i < numberRows_; i++) { char dropped = static_cast(rowsDropped[i]); rowsDropped_[i] = dropped; if (dropped == 2) { /*dropped this time*/ rowsDropped[newDropped++] = i; rowsDropped_[i] = 0; } } numberRowsDropped_ = newDropped; newDropped = -(2 + newDropped); } } else { /* KKT*/ CoinPackedMatrix * quadratic = NULL; ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(model_->objectiveAsObject())); if (quadraticObj) quadratic = quadraticObj->quadraticObjective(); /* matrix*/ int numberRowsModel = model_->numberRows(); int numberColumns = model_->numberColumns(); int numberTotal = numberColumns + numberRowsModel; longDouble * work = sparseFactor_; work--; /* skip diagonal*/ int addOffset = numberRows_ - 1; int iColumn; if (!quadratic) { for (iColumn = 0; iColumn < numberColumns; iColumn++) { CoinWorkDouble value = diagonal[iColumn]; if (CoinAbs(value) > 1.0e-100) { value = 1.0 / value; largest = CoinMax(largest, CoinAbs(value)); diagonal_[iColumn] = -value; CoinBigIndex start = columnStart[iColumn]; CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn]; for (CoinBigIndex j = start; j < end; j++) { /*choleskyRow_[numberElements]=row[j]+numberTotal;*/ /*sparseFactor_[numberElements++]=element[j];*/ work[row[j] + numberTotal] = element[j]; largest = CoinMax(largest, CoinAbs(element[j])); } } else { diagonal_[iColumn] = -value; } addOffset--; work += addOffset; } } else { /* Quadratic*/ const int * columnQuadratic = quadratic->getIndices(); const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts(); const int * columnQuadraticLength = quadratic->getVectorLengths(); const double * quadraticElement = quadratic->getElements(); for (iColumn = 0; iColumn < numberColumns; iColumn++) { CoinWorkDouble value = diagonal[iColumn]; CoinBigIndex j; if (CoinAbs(value) > 1.0e-100) { value = 1.0 / value; for (j = columnQuadraticStart[iColumn]; j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) { int jColumn = columnQuadratic[j]; if (jColumn > iColumn) { work[jColumn] = -quadraticElement[j]; } else if (iColumn == jColumn) { value += quadraticElement[j]; } } largest = CoinMax(largest, CoinAbs(value)); diagonal_[iColumn] = -value; CoinBigIndex start = columnStart[iColumn]; CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn]; for (j = start; j < end; j++) { work[row[j] + numberTotal] = element[j]; largest = CoinMax(largest, CoinAbs(element[j])); } } else { value = 1.0e100; diagonal_[iColumn] = -value; } addOffset--; work += addOffset; } } /* slacks*/ for (iColumn = numberColumns; iColumn < numberTotal; iColumn++) { CoinWorkDouble value = diagonal[iColumn]; if (CoinAbs(value) > 1.0e-100) { value = 1.0 / value; largest = CoinMax(largest, CoinAbs(value)); } else { value = 1.0e100; } diagonal_[iColumn] = -value; work[iColumn-numberColumns+numberTotal] = -1.0; addOffset--; work += addOffset; } /* Finish diagonal*/ for (iRow = 0; iRow < numberRowsModel; iRow++) { diagonal_[iRow+numberTotal] = delta2; } /*check sizes*/ largest *= 1.0e-20; largest = CoinMin(largest, CHOL_SMALL_VALUE); doubleParameters_[10] = CoinMax(1.0e-20, largest); integerParameters_[20] = 0; doubleParameters_[3] = 0.0; doubleParameters_[4] = COIN_DBL_MAX; /* Set up LDL cutoff*/ integerParameters_[34] = numberTotal; /* KKT*/ int * rowsDropped2 = new int[numberRows_]; CoinZeroN(rowsDropped2, numberRows_); #ifdef CHOL_COMPARE if (numberRows_ < 200) factorizePart3(rowsDropped2); #endif factorizePart2(rowsDropped2); newDropped = integerParameters_[20]; largest = doubleParameters_[3]; smallest = doubleParameters_[4]; if (model_->messageHandler()->logLevel() > 1) COIN_DETAIL_PRINT(std::cout << "Cholesky - largest " << largest << " smallest " << smallest << std::endl); choleskyCondition_ = largest / smallest; /* Should save adjustments in ..R_*/ int n1 = 0, n2 = 0; CoinWorkDouble * primalR = model_->primalR(); CoinWorkDouble * dualR = model_->dualR(); for (iRow = 0; iRow < numberTotal; iRow++) { if (rowsDropped2[iRow]) { n1++; /*printf("row region1 %d dropped\n",iRow);*/ /*rowsDropped_[iRow]=1;*/ rowsDropped_[iRow] = 0; primalR[iRow] = doubleParameters_[20]; } else { rowsDropped_[iRow] = 0; primalR[iRow] = 0.0; } } for (; iRow < numberRows_; iRow++) { if (rowsDropped2[iRow]) { n2++; /*printf("row region2 %d dropped\n",iRow);*/ /*rowsDropped_[iRow]=1;*/ rowsDropped_[iRow] = 0; dualR[iRow-numberTotal] = doubleParameters_[34]; } else { rowsDropped_[iRow] = 0; dualR[iRow-numberTotal] = 0.0; } } } return 0; } /* Factorize - filling in rowsDropped and returning number dropped */ void ClpCholeskyDense::factorizePart3( int * rowsDropped) { int iColumn; longDouble * xx = sparseFactor_; longDouble * yy = diagonal_; diagonal_ = sparseFactor_ + 40000; sparseFactor_ = diagonal_ + numberRows_; /*memcpy(sparseFactor_,xx,sizeFactor_*sizeof(double));*/ CoinMemcpyN(xx, 40000, sparseFactor_); CoinMemcpyN(yy, numberRows_, diagonal_); int numberDropped = 0; CoinWorkDouble largest = 0.0; CoinWorkDouble smallest = COIN_DBL_MAX; double dropValue = doubleParameters_[10]; int firstPositive = integerParameters_[34]; longDouble * work = sparseFactor_; /* Allow for triangular*/ int addOffset = numberRows_ - 1; work--; for (iColumn = 0; iColumn < numberRows_; iColumn++) { int iRow; int addOffsetNow = numberRows_ - 1;; longDouble * workNow = sparseFactor_ - 1 + iColumn; CoinWorkDouble diagonalValue = diagonal_[iColumn]; for (iRow = 0; iRow < iColumn; iRow++) { double aj = *workNow; addOffsetNow--; workNow += addOffsetNow; diagonalValue -= aj * aj * workDouble_[iRow]; } bool dropColumn = false; if (iColumn < firstPositive) { /* must be negative*/ if (diagonalValue <= -dropValue) { smallest = CoinMin(smallest, -diagonalValue); largest = CoinMax(largest, -diagonalValue); workDouble_[iColumn] = diagonalValue; diagonalValue = 1.0 / diagonalValue; } else { dropColumn = true; workDouble_[iColumn] = -1.0e100; diagonalValue = 0.0; integerParameters_[20]++; } } else { /* must be positive*/ if (diagonalValue >= dropValue) { smallest = CoinMin(smallest, diagonalValue); largest = CoinMax(largest, diagonalValue); workDouble_[iColumn] = diagonalValue; diagonalValue = 1.0 / diagonalValue; } else { dropColumn = true; workDouble_[iColumn] = 1.0e100; diagonalValue = 0.0; integerParameters_[20]++; } } if (!dropColumn) { diagonal_[iColumn] = diagonalValue; for (iRow = iColumn + 1; iRow < numberRows_; iRow++) { double value = work[iRow]; workNow = sparseFactor_ - 1; int addOffsetNow = numberRows_ - 1;; for (int jColumn = 0; jColumn < iColumn; jColumn++) { double aj = workNow[iColumn]; double multiplier = workDouble_[jColumn]; double ai = workNow[iRow]; addOffsetNow--; workNow += addOffsetNow; value -= aj * ai * multiplier; } work[iRow] = value * diagonalValue; } } else { /* drop column*/ rowsDropped[iColumn] = 2; numberDropped++; diagonal_[iColumn] = 0.0; for (iRow = iColumn + 1; iRow < numberRows_; iRow++) { work[iRow] = 0.0; } } addOffset--; work += addOffset; } doubleParameters_[3] = largest; doubleParameters_[4] = smallest; integerParameters_[20] = numberDropped; sparseFactor_ = xx; diagonal_ = yy; } /*#define POS_DEBUG*/ #ifdef POS_DEBUG static int counter = 0; int ClpCholeskyDense::bNumber(const longDouble * array, int &iRow, int &iCol) { int numberBlocks = (numberRows_ + BLOCK - 1) >> BLOCKSHIFT; longDouble * a = sparseFactor_ + BLOCKSQ * numberBlocks; int k = array - a; assert ((k % BLOCKSQ) == 0); iCol = 0; int kk = k >> BLOCKSQSHIFT; /*printf("%d %d %d %d\n",k,kk,BLOCKSQ,BLOCKSQSHIFT);*/ k = kk; while (k >= numberBlocks) { iCol++; k -= numberBlocks; numberBlocks--; } iRow = iCol + k; counter++; if(counter > 100000) exit(77); return kk; } #endif /* Factorize - filling in rowsDropped and returning number dropped */ void ClpCholeskyDense::factorizePart2( int * rowsDropped) { int iColumn; /*longDouble * xx = sparseFactor_;*/ /*longDouble * yy = diagonal_;*/ /*diagonal_ = sparseFactor_ + 40000;*/ /*sparseFactor_=diagonal_ + numberRows_;*/ /*memcpy(sparseFactor_,xx,sizeFactor_*sizeof(double));*/ /*memcpy(sparseFactor_,xx,40000*sizeof(double));*/ /*memcpy(diagonal_,yy,numberRows_*sizeof(double));*/ int numberBlocks = (numberRows_ + BLOCK - 1) >> BLOCKSHIFT; /* later align on boundary*/ longDouble * a = sparseFactor_ + BLOCKSQ * numberBlocks; int n = numberRows_; int nRound = numberRows_ & (~(BLOCK - 1)); /* adjust if exact*/ if (nRound == n) nRound -= BLOCK; int sizeLastBlock = n - nRound; int get = n * (n - 1) / 2; /* ? as no diagonal*/ int block = numberBlocks * (numberBlocks + 1) / 2; int ifOdd; int rowLast; if (sizeLastBlock != BLOCK) { longDouble * aa = &a[(block-1)*BLOCKSQ]; rowLast = nRound - 1; ifOdd = 1; int put = BLOCKSQ; /* do last separately*/ put -= (BLOCK - sizeLastBlock) * (BLOCK + 1); for (iColumn = numberRows_ - 1; iColumn >= nRound; iColumn--) { int put2 = put; put -= BLOCK; for (int iRow = numberRows_ - 1; iRow > iColumn; iRow--) { aa[--put2] = sparseFactor_[--get]; assert (aa + put2 >= sparseFactor_ + get); } /* save diagonal as well*/ aa[--put2] = diagonal_[iColumn]; } n = nRound; block--; } else { /* exact fit*/ rowLast = numberRows_ - 1; ifOdd = 0; } /* Now main loop*/ int nBlock = 0; for (; n > 0; n -= BLOCK) { longDouble * aa = &a[(block-1)*BLOCKSQ]; longDouble * aaLast = NULL; int put = BLOCKSQ; int putLast = 0; /* see if we have small block*/ if (ifOdd) { aaLast = &a[(block-1)*BLOCKSQ]; aa = aaLast - BLOCKSQ; putLast = BLOCKSQ - BLOCK + sizeLastBlock; } for (iColumn = n - 1; iColumn >= n - BLOCK; iColumn--) { if (aaLast) { /* last bit*/ for (int iRow = numberRows_ - 1; iRow > rowLast; iRow--) { aaLast[--putLast] = sparseFactor_[--get]; assert (aaLast + putLast >= sparseFactor_ + get); } putLast -= BLOCK - sizeLastBlock; } longDouble * aPut = aa; int j = rowLast; for (int jBlock = 0; jBlock <= nBlock; jBlock++) { int put2 = put; int last = CoinMax(j - BLOCK, iColumn); for (int iRow = j; iRow > last; iRow--) { aPut[--put2] = sparseFactor_[--get]; assert (aPut + put2 >= sparseFactor_ + get); } if (j - BLOCK < iColumn) { /* save diagonal as well*/ aPut[--put2] = diagonal_[iColumn]; } j -= BLOCK; aPut -= BLOCKSQ; } put -= BLOCK; } nBlock++; block -= nBlock + ifOdd; } ClpCholeskyDenseC info; info.diagonal_ = diagonal_; info.doubleParameters_[0] = doubleParameters_[10]; info.integerParameters_[0] = integerParameters_[34]; #ifndef CLP_CILK ClpCholeskyCfactor(&info, a, numberRows_, numberBlocks, diagonal_, workDouble_, rowsDropped); #else info.a = a; info.n = numberRows_; info.numberBlocks = numberBlocks; info.work = workDouble_; info.rowsDropped = rowsDropped; info.integerParameters_[1] = model_->numberThreads(); ClpCholeskySpawn(&info); #endif double largest = 0.0; double smallest = COIN_DBL_MAX; int numberDropped = 0; for (int i = 0; i < numberRows_; i++) { if (diagonal_[i]) { largest = CoinMax(largest, CoinAbs(diagonal_[i])); smallest = CoinMin(smallest, CoinAbs(diagonal_[i])); } else { numberDropped++; } } doubleParameters_[3] = CoinMax(doubleParameters_[3], 1.0 / smallest); doubleParameters_[4] = CoinMin(doubleParameters_[4], 1.0 / largest); integerParameters_[20] += numberDropped; } /* Non leaf recursive factor*/ void ClpCholeskyCfactor(ClpCholeskyDenseC * thisStruct, longDouble * a, int n, int numberBlocks, longDouble * diagonal, longDouble * work, int * rowsDropped) { if (n <= BLOCK) { ClpCholeskyCfactorLeaf(thisStruct, a, n, diagonal, work, rowsDropped); } else { int nb = number_blocks((n + 1) >> 1); int nThis = number_rows(nb); longDouble * aother; int nLeft = n - nThis; int nintri = (nb * (nb + 1)) >> 1; int nbelow = (numberBlocks - nb) * nb; ClpCholeskyCfactor(thisStruct, a, nThis, numberBlocks, diagonal, work, rowsDropped); ClpCholeskyCtriRec(thisStruct, a, nThis, a + number_entries(nb), diagonal, work, nLeft, nb, 0, numberBlocks); aother = a + number_entries(nintri + nbelow); ClpCholeskyCrecTri(thisStruct, a + number_entries(nb), nLeft, nThis, nb, 0, aother, diagonal, work, numberBlocks); ClpCholeskyCfactor(thisStruct, aother, nLeft, numberBlocks - nb, diagonal + nThis, work + nThis, rowsDropped); } } /* Non leaf recursive triangle rectangle update*/ void ClpCholeskyCtriRec(ClpCholeskyDenseC * thisStruct, longDouble * aTri, int nThis, longDouble * aUnder, longDouble * diagonal, longDouble * work, int nLeft, int iBlock, int jBlock, int numberBlocks) { if (nThis <= BLOCK && nLeft <= BLOCK) { ClpCholeskyCtriRecLeaf(/*thisStruct,*/ aTri, aUnder, diagonal, work, nLeft); } else if (nThis < nLeft) { int nb = number_blocks((nLeft + 1) >> 1); int nLeft2 = number_rows(nb); cilk_spawn ClpCholeskyCtriRec(thisStruct, aTri, nThis, aUnder, diagonal, work, nLeft2, iBlock, jBlock, numberBlocks); ClpCholeskyCtriRec(thisStruct, aTri, nThis, aUnder + number_entries(nb), diagonal, work, nLeft - nLeft2, iBlock + nb, jBlock, numberBlocks); cilk_sync; } else { int nb = number_blocks((nThis + 1) >> 1); int nThis2 = number_rows(nb); longDouble * aother; int kBlock = jBlock + nb; int i; int nintri = (nb * (nb + 1)) >> 1; int nbelow = (numberBlocks - nb) * nb; ClpCholeskyCtriRec(thisStruct, aTri, nThis2, aUnder, diagonal, work, nLeft, iBlock, jBlock, numberBlocks); /* and rectangular update */ i = ((numberBlocks - jBlock) * (numberBlocks - jBlock - 1) - (numberBlocks - jBlock - nb) * (numberBlocks - jBlock - nb - 1)) >> 1; aother = aUnder + number_entries(i); ClpCholeskyCrecRec(thisStruct, aTri + number_entries(nb), nThis - nThis2, nLeft, nThis2, aUnder, aother, work, kBlock, jBlock, numberBlocks); ClpCholeskyCtriRec(thisStruct, aTri + number_entries(nintri + nbelow), nThis - nThis2, aother, diagonal + nThis2, work + nThis2, nLeft, iBlock - nb, kBlock - nb, numberBlocks - nb); } } /* Non leaf recursive rectangle triangle update*/ void ClpCholeskyCrecTri(ClpCholeskyDenseC * thisStruct, longDouble * aUnder, int nTri, int nDo, int iBlock, int jBlock, longDouble * aTri, longDouble * diagonal, longDouble * work, int numberBlocks) { if (nTri <= BLOCK && nDo <= BLOCK) { ClpCholeskyCrecTriLeaf(/*thisStruct,*/ aUnder, aTri,/*diagonal,*/work, nTri); } else if (nTri < nDo) { int nb = number_blocks((nDo + 1) >> 1); int nDo2 = number_rows(nb); longDouble * aother; int i; ClpCholeskyCrecTri(thisStruct, aUnder, nTri, nDo2, iBlock, jBlock, aTri, diagonal, work, numberBlocks); i = ((numberBlocks - jBlock) * (numberBlocks - jBlock - 1) - (numberBlocks - jBlock - nb) * (numberBlocks - jBlock - nb - 1)) >> 1; aother = aUnder + number_entries(i); ClpCholeskyCrecTri(thisStruct, aother, nTri, nDo - nDo2, iBlock - nb, jBlock, aTri, diagonal + nDo2, work + nDo2, numberBlocks - nb); } else { int nb = number_blocks((nTri + 1) >> 1); int nTri2 = number_rows(nb); longDouble * aother; int i; cilk_spawn ClpCholeskyCrecTri(thisStruct, aUnder, nTri2, nDo, iBlock, jBlock, aTri, diagonal, work, numberBlocks); /* and rectangular update */ i = ((numberBlocks - iBlock) * (numberBlocks - iBlock + 1) - (numberBlocks - iBlock - nb) * (numberBlocks - iBlock - nb + 1)) >> 1; aother = aTri + number_entries(nb); ClpCholeskyCrecRec(thisStruct, aUnder, nTri2, nTri - nTri2, nDo, aUnder + number_entries(nb), aother, work, iBlock, jBlock, numberBlocks); ClpCholeskyCrecTri(thisStruct, aUnder + number_entries(nb), nTri - nTri2, nDo, iBlock + nb, jBlock, aTri + number_entries(i), diagonal, work, numberBlocks); cilk_sync; } } /* Non leaf recursive rectangle rectangle update, nUnder is number of rows in iBlock, nUnderK is number of rows in kBlock */ void ClpCholeskyCrecRec(ClpCholeskyDenseC * thisStruct, longDouble * above, int nUnder, int nUnderK, int nDo, longDouble * aUnder, longDouble *aOther, longDouble * work, int iBlock, int jBlock, int numberBlocks) { if (nDo <= BLOCK && nUnder <= BLOCK && nUnderK <= BLOCK) { assert (nDo == BLOCK && nUnder == BLOCK); ClpCholeskyCrecRecLeaf(/*thisStruct,*/ above , aUnder , aOther, work, nUnderK); } else if (nDo <= nUnderK && nUnder <= nUnderK) { int nb = number_blocks((nUnderK + 1) >> 1); int nUnder2 = number_rows(nb); cilk_spawn ClpCholeskyCrecRec(thisStruct, above, nUnder, nUnder2, nDo, aUnder, aOther, work, iBlock, jBlock, numberBlocks); ClpCholeskyCrecRec(thisStruct, above, nUnder, nUnderK - nUnder2, nDo, aUnder + number_entries(nb), aOther + number_entries(nb), work, iBlock, jBlock, numberBlocks); cilk_sync; } else if (nUnderK <= nDo && nUnder <= nDo) { int nb = number_blocks((nDo + 1) >> 1); int nDo2 = number_rows(nb); int i; ClpCholeskyCrecRec(thisStruct, above, nUnder, nUnderK, nDo2, aUnder, aOther, work, iBlock, jBlock, numberBlocks); i = ((numberBlocks - jBlock) * (numberBlocks - jBlock - 1) - (numberBlocks - jBlock - nb) * (numberBlocks - jBlock - nb - 1)) >> 1; ClpCholeskyCrecRec(thisStruct, above + number_entries(i), nUnder, nUnderK, nDo - nDo2, aUnder + number_entries(i), aOther, work + nDo2, iBlock - nb, jBlock, numberBlocks - nb); } else { int nb = number_blocks((nUnder + 1) >> 1); int nUnder2 = number_rows(nb); int i; cilk_spawn ClpCholeskyCrecRec(thisStruct, above, nUnder2, nUnderK, nDo, aUnder, aOther, work, iBlock, jBlock, numberBlocks); i = ((numberBlocks - iBlock) * (numberBlocks - iBlock - 1) - (numberBlocks - iBlock - nb) * (numberBlocks - iBlock - nb - 1)) >> 1; ClpCholeskyCrecRec(thisStruct, above + number_entries(nb), nUnder - nUnder2, nUnderK, nDo, aUnder, aOther + number_entries(i), work, iBlock + nb, jBlock, numberBlocks); cilk_sync; } } /* Leaf recursive factor*/ void ClpCholeskyCfactorLeaf(ClpCholeskyDenseC * thisStruct, longDouble * a, int n, longDouble * diagonal, longDouble * work, int * rowsDropped) { double dropValue = thisStruct->doubleParameters_[0]; int firstPositive = thisStruct->integerParameters_[0]; int rowOffset = static_cast(diagonal - thisStruct->diagonal_); int i, j, k; CoinWorkDouble t00, temp1; longDouble * aa; aa = a - BLOCK; for (j = 0; j < n; j ++) { bool dropColumn; CoinWorkDouble useT00; aa += BLOCK; t00 = aa[j]; for (k = 0; k < j; ++k) { CoinWorkDouble multiplier = work[k]; t00 -= a[j + k * BLOCK] * a[j + k * BLOCK] * multiplier; } dropColumn = false; useT00 = t00; if (j + rowOffset < firstPositive) { /* must be negative*/ if (t00 <= -dropValue) { /*aa[j]=t00;*/ t00 = 1.0 / t00; } else { dropColumn = true; /*aa[j]=-1.0e100;*/ useT00 = -1.0e-100; t00 = 0.0; } } else { /* must be positive*/ if (t00 >= dropValue) { /*aa[j]=t00;*/ t00 = 1.0 / t00; } else { dropColumn = true; /*aa[j]=1.0e100;*/ useT00 = 1.0e-100; t00 = 0.0; } } if (!dropColumn) { diagonal[j] = t00; work[j] = useT00; temp1 = t00; for (i = j + 1; i < n; i++) { t00 = aa[i]; for (k = 0; k < j; ++k) { CoinWorkDouble multiplier = work[k]; t00 -= a[i + k * BLOCK] * a[j + k * BLOCK] * multiplier; } aa[i] = t00 * temp1; } } else { /* drop column*/ rowsDropped[j+rowOffset] = 2; diagonal[j] = 0.0; /*aa[j]=1.0e100;*/ work[j] = 1.0e100; for (i = j + 1; i < n; i++) { aa[i] = 0.0; } } } } /* Leaf recursive triangle rectangle update*/ void ClpCholeskyCtriRecLeaf(/*ClpCholeskyDenseC * thisStruct,*/ longDouble * aTri, longDouble * aUnder, longDouble * diagonal, longDouble * work, int nUnder) { #ifdef POS_DEBUG int iru, icu; int iu = bNumber(aUnder, iru, icu); int irt, ict; int it = bNumber(aTri, irt, ict); /*printf("%d %d\n",iu,it);*/ printf("trirecleaf under (%d,%d), tri (%d,%d)\n", iru, icu, irt, ict); assert (diagonal == thisStruct->diagonal_ + ict * BLOCK); #endif int j; longDouble * aa; #ifdef BLOCKUNROLL if (nUnder == BLOCK) { aa = aTri - 2 * BLOCK; for (j = 0; j < BLOCK; j += 2) { int i; CoinWorkDouble temp0 = diagonal[j]; CoinWorkDouble temp1 = diagonal[j+1]; aa += 2 * BLOCK; for ( i = 0; i < BLOCK; i += 2) { CoinWorkDouble at1; CoinWorkDouble t00 = aUnder[i+j*BLOCK]; CoinWorkDouble t10 = aUnder[i+ BLOCK + j*BLOCK]; CoinWorkDouble t01 = aUnder[i+1+j*BLOCK]; CoinWorkDouble t11 = aUnder[i+1+ BLOCK + j*BLOCK]; int k; for (k = 0; k < j; ++k) { CoinWorkDouble multiplier = work[k]; CoinWorkDouble au0 = aUnder[i + k * BLOCK] * multiplier; CoinWorkDouble au1 = aUnder[i + 1 + k * BLOCK] * multiplier; CoinWorkDouble at0 = aTri[j + k * BLOCK]; CoinWorkDouble at1 = aTri[j + 1 + k * BLOCK]; t00 -= au0 * at0; t10 -= au0 * at1; t01 -= au1 * at0; t11 -= au1 * at1; } t00 *= temp0; at1 = aTri[j + 1 + j*BLOCK] * work[j]; t10 -= t00 * at1; t01 *= temp0; t11 -= t01 * at1; aUnder[i+j*BLOCK] = t00; aUnder[i+1+j*BLOCK] = t01; aUnder[i+ BLOCK + j*BLOCK] = t10 * temp1; aUnder[i+1+ BLOCK + j*BLOCK] = t11 * temp1; } } } else { #endif aa = aTri - BLOCK; for (j = 0; j < BLOCK; j ++) { int i; CoinWorkDouble temp1 = diagonal[j]; aa += BLOCK; for (i = 0; i < nUnder; i++) { int k; CoinWorkDouble t00 = aUnder[i+j*BLOCK]; for ( k = 0; k < j; ++k) { CoinWorkDouble multiplier = work[k]; t00 -= aUnder[i + k * BLOCK] * aTri[j + k * BLOCK] * multiplier; } aUnder[i+j*BLOCK] = t00 * temp1; } } #ifdef BLOCKUNROLL } #endif } /* Leaf recursive rectangle triangle update*/ void ClpCholeskyCrecTriLeaf(/*ClpCholeskyDenseC * thisStruct,*/ longDouble * aUnder, longDouble * aTri, /*longDouble * diagonal,*/ longDouble * work, int nUnder) { #ifdef POS_DEBUG int iru, icu; int iu = bNumber(aUnder, iru, icu); int irt, ict; int it = bNumber(aTri, irt, ict); /*printf("%d %d\n",iu,it);*/ printf("rectrileaf under (%d,%d), tri (%d,%d)\n", iru, icu, irt, ict); assert (diagonal == thisStruct->diagonal_ + icu * BLOCK); #endif int i, j, k; CoinWorkDouble t00; longDouble * aa; #ifdef BLOCKUNROLL if (nUnder == BLOCK) { longDouble * aUnder2 = aUnder - 2; aa = aTri - 2 * BLOCK; for (j = 0; j < BLOCK; j += 2) { CoinWorkDouble t00, t01, t10, t11; aa += 2 * BLOCK; aUnder2 += 2; t00 = aa[j]; t01 = aa[j+1]; t10 = aa[j+1+BLOCK]; for (k = 0; k < BLOCK; ++k) { CoinWorkDouble multiplier = work[k]; CoinWorkDouble a0 = aUnder2[k * BLOCK]; CoinWorkDouble a1 = aUnder2[1 + k * BLOCK]; CoinWorkDouble x0 = a0 * multiplier; CoinWorkDouble x1 = a1 * multiplier; t00 -= a0 * x0; t01 -= a1 * x0; t10 -= a1 * x1; } aa[j] = t00; aa[j+1] = t01; aa[j+1+BLOCK] = t10; for (i = j + 2; i < BLOCK; i += 2) { t00 = aa[i]; t01 = aa[i+BLOCK]; t10 = aa[i+1]; t11 = aa[i+1+BLOCK]; for (k = 0; k < BLOCK; ++k) { CoinWorkDouble multiplier = work[k]; CoinWorkDouble a0 = aUnder2[k * BLOCK] * multiplier; CoinWorkDouble a1 = aUnder2[1 + k * BLOCK] * multiplier; t00 -= aUnder[i + k * BLOCK] * a0; t01 -= aUnder[i + k * BLOCK] * a1; t10 -= aUnder[i + 1 + k * BLOCK] * a0; t11 -= aUnder[i + 1 + k * BLOCK] * a1; } aa[i] = t00; aa[i+BLOCK] = t01; aa[i+1] = t10; aa[i+1+BLOCK] = t11; } } } else { #endif aa = aTri - BLOCK; for (j = 0; j < nUnder; j ++) { aa += BLOCK; for (i = j; i < nUnder; i++) { t00 = aa[i]; for (k = 0; k < BLOCK; ++k) { CoinWorkDouble multiplier = work[k]; t00 -= aUnder[i + k * BLOCK] * aUnder[j + k * BLOCK] * multiplier; } aa[i] = t00; } } #ifdef BLOCKUNROLL } #endif } /* Leaf recursive rectangle rectangle update, nUnder is number of rows in iBlock, nUnderK is number of rows in kBlock */ void ClpCholeskyCrecRecLeaf(/*ClpCholeskyDenseC * thisStruct,*/ const longDouble * COIN_RESTRICT above, const longDouble * COIN_RESTRICT aUnder, longDouble * COIN_RESTRICT aOther, const longDouble * COIN_RESTRICT work, int nUnder) { #ifdef POS_DEBUG int ira, ica; int ia = bNumber(above, ira, ica); int iru, icu; int iu = bNumber(aUnder, iru, icu); int iro, ico; int io = bNumber(aOther, iro, ico); /*printf("%d %d %d\n",ia,iu,io);*/ printf("recrecleaf above (%d,%d), under (%d,%d), other (%d,%d)\n", ira, ica, iru, icu, iro, ico); #endif int i, j, k; longDouble * aa; #ifdef BLOCKUNROLL aa = aOther - 4 * BLOCK; if (nUnder == BLOCK) { /*#define INTEL*/ #ifdef INTEL aa += 2 * BLOCK; for (j = 0; j < BLOCK; j += 2) { aa += 2 * BLOCK; for (i = 0; i < BLOCK; i += 2) { CoinWorkDouble t00 = aa[i+0*BLOCK]; CoinWorkDouble t10 = aa[i+1*BLOCK]; CoinWorkDouble t01 = aa[i+1+0*BLOCK]; CoinWorkDouble t11 = aa[i+1+1*BLOCK]; for (k = 0; k < BLOCK; k++) { CoinWorkDouble multiplier = work[k]; CoinWorkDouble a00 = aUnder[i+k*BLOCK] * multiplier; CoinWorkDouble a01 = aUnder[i+1+k*BLOCK] * multiplier; t00 -= a00 * above[j + 0 + k * BLOCK]; t10 -= a00 * above[j + 1 + k * BLOCK]; t01 -= a01 * above[j + 0 + k * BLOCK]; t11 -= a01 * above[j + 1 + k * BLOCK]; } aa[i+0*BLOCK] = t00; aa[i+1*BLOCK] = t10; aa[i+1+0*BLOCK] = t01; aa[i+1+1*BLOCK] = t11; } } #else for (j = 0; j < BLOCK; j += 4) { aa += 4 * BLOCK; for (i = 0; i < BLOCK; i += 4) { CoinWorkDouble t00 = aa[i+0+0*BLOCK]; CoinWorkDouble t10 = aa[i+0+1*BLOCK]; CoinWorkDouble t20 = aa[i+0+2*BLOCK]; CoinWorkDouble t30 = aa[i+0+3*BLOCK]; CoinWorkDouble t01 = aa[i+1+0*BLOCK]; CoinWorkDouble t11 = aa[i+1+1*BLOCK]; CoinWorkDouble t21 = aa[i+1+2*BLOCK]; CoinWorkDouble t31 = aa[i+1+3*BLOCK]; CoinWorkDouble t02 = aa[i+2+0*BLOCK]; CoinWorkDouble t12 = aa[i+2+1*BLOCK]; CoinWorkDouble t22 = aa[i+2+2*BLOCK]; CoinWorkDouble t32 = aa[i+2+3*BLOCK]; CoinWorkDouble t03 = aa[i+3+0*BLOCK]; CoinWorkDouble t13 = aa[i+3+1*BLOCK]; CoinWorkDouble t23 = aa[i+3+2*BLOCK]; CoinWorkDouble t33 = aa[i+3+3*BLOCK]; const longDouble * COIN_RESTRICT aUnderNow = aUnder + i; const longDouble * COIN_RESTRICT aboveNow = above + j; for (k = 0; k < BLOCK; k++) { CoinWorkDouble multiplier = work[k]; CoinWorkDouble a00 = aUnderNow[0] * multiplier; CoinWorkDouble a01 = aUnderNow[1] * multiplier; CoinWorkDouble a02 = aUnderNow[2] * multiplier; CoinWorkDouble a03 = aUnderNow[3] * multiplier; t00 -= a00 * aboveNow[0]; t10 -= a00 * aboveNow[1]; t20 -= a00 * aboveNow[2]; t30 -= a00 * aboveNow[3]; t01 -= a01 * aboveNow[0]; t11 -= a01 * aboveNow[1]; t21 -= a01 * aboveNow[2]; t31 -= a01 * aboveNow[3]; t02 -= a02 * aboveNow[0]; t12 -= a02 * aboveNow[1]; t22 -= a02 * aboveNow[2]; t32 -= a02 * aboveNow[3]; t03 -= a03 * aboveNow[0]; t13 -= a03 * aboveNow[1]; t23 -= a03 * aboveNow[2]; t33 -= a03 * aboveNow[3]; aUnderNow += BLOCK; aboveNow += BLOCK; } aa[i+0+0*BLOCK] = t00; aa[i+0+1*BLOCK] = t10; aa[i+0+2*BLOCK] = t20; aa[i+0+3*BLOCK] = t30; aa[i+1+0*BLOCK] = t01; aa[i+1+1*BLOCK] = t11; aa[i+1+2*BLOCK] = t21; aa[i+1+3*BLOCK] = t31; aa[i+2+0*BLOCK] = t02; aa[i+2+1*BLOCK] = t12; aa[i+2+2*BLOCK] = t22; aa[i+2+3*BLOCK] = t32; aa[i+3+0*BLOCK] = t03; aa[i+3+1*BLOCK] = t13; aa[i+3+2*BLOCK] = t23; aa[i+3+3*BLOCK] = t33; } } #endif } else { int odd = nUnder & 1; int n = nUnder - odd; for (j = 0; j < BLOCK; j += 4) { aa += 4 * BLOCK; for (i = 0; i < n; i += 2) { CoinWorkDouble t00 = aa[i+0*BLOCK]; CoinWorkDouble t10 = aa[i+1*BLOCK]; CoinWorkDouble t20 = aa[i+2*BLOCK]; CoinWorkDouble t30 = aa[i+3*BLOCK]; CoinWorkDouble t01 = aa[i+1+0*BLOCK]; CoinWorkDouble t11 = aa[i+1+1*BLOCK]; CoinWorkDouble t21 = aa[i+1+2*BLOCK]; CoinWorkDouble t31 = aa[i+1+3*BLOCK]; const longDouble * COIN_RESTRICT aUnderNow = aUnder + i; const longDouble * COIN_RESTRICT aboveNow = above + j; for (k = 0; k < BLOCK; k++) { CoinWorkDouble multiplier = work[k]; CoinWorkDouble a00 = aUnderNow[0] * multiplier; CoinWorkDouble a01 = aUnderNow[1] * multiplier; t00 -= a00 * aboveNow[0]; t10 -= a00 * aboveNow[1]; t20 -= a00 * aboveNow[2]; t30 -= a00 * aboveNow[3]; t01 -= a01 * aboveNow[0]; t11 -= a01 * aboveNow[1]; t21 -= a01 * aboveNow[2]; t31 -= a01 * aboveNow[3]; aUnderNow += BLOCK; aboveNow += BLOCK; } aa[i+0*BLOCK] = t00; aa[i+1*BLOCK] = t10; aa[i+2*BLOCK] = t20; aa[i+3*BLOCK] = t30; aa[i+1+0*BLOCK] = t01; aa[i+1+1*BLOCK] = t11; aa[i+1+2*BLOCK] = t21; aa[i+1+3*BLOCK] = t31; } if (odd) { CoinWorkDouble t0 = aa[n+0*BLOCK]; CoinWorkDouble t1 = aa[n+1*BLOCK]; CoinWorkDouble t2 = aa[n+2*BLOCK]; CoinWorkDouble t3 = aa[n+3*BLOCK]; CoinWorkDouble a0; for (k = 0; k < BLOCK; k++) { a0 = aUnder[n+k*BLOCK] * work[k]; t0 -= a0 * above[j + 0 + k * BLOCK]; t1 -= a0 * above[j + 1 + k * BLOCK]; t2 -= a0 * above[j + 2 + k * BLOCK]; t3 -= a0 * above[j + 3 + k * BLOCK]; } aa[n+0*BLOCK] = t0; aa[n+1*BLOCK] = t1; aa[n+2*BLOCK] = t2; aa[n+3*BLOCK] = t3; } } } #else aa = aOther - BLOCK; for (j = 0; j < BLOCK; j ++) { aa += BLOCK; for (i = 0; i < nUnder; i++) { CoinWorkDouble t00 = aa[i+0*BLOCK]; for (k = 0; k < BLOCK; k++) { CoinWorkDouble a00 = aUnder[i+k*BLOCK] * work[k]; t00 -= a00 * above[j + k * BLOCK]; } aa[i] = t00; } } #endif } /* Uses factorization to solve. */ void ClpCholeskyDense::solve (CoinWorkDouble * region) { #ifdef CHOL_COMPARE double * region2 = NULL; if (numberRows_ < 200) { longDouble * xx = sparseFactor_; longDouble * yy = diagonal_; diagonal_ = sparseFactor_ + 40000; sparseFactor_ = diagonal_ + numberRows_; region2 = ClpCopyOfArray(region, numberRows_); int iRow, iColumn; int addOffset = numberRows_ - 1; longDouble * work = sparseFactor_ - 1; for (iColumn = 0; iColumn < numberRows_; iColumn++) { double value = region2[iColumn]; for (iRow = iColumn + 1; iRow < numberRows_; iRow++) region2[iRow] -= value * work[iRow]; addOffset--; work += addOffset; } for (iColumn = numberRows_ - 1; iColumn >= 0; iColumn--) { double value = region2[iColumn] * diagonal_[iColumn]; work -= addOffset; addOffset++; for (iRow = iColumn + 1; iRow < numberRows_; iRow++) value -= region2[iRow] * work[iRow]; region2[iColumn] = value; } sparseFactor_ = xx; diagonal_ = yy; } #endif /*longDouble * xx = sparseFactor_;*/ /*longDouble * yy = diagonal_;*/ /*diagonal_ = sparseFactor_ + 40000;*/ /*sparseFactor_=diagonal_ + numberRows_;*/ int numberBlocks = (numberRows_ + BLOCK - 1) >> BLOCKSHIFT; /* later align on boundary*/ longDouble * a = sparseFactor_ + BLOCKSQ * numberBlocks; int iBlock; longDouble * aa = a; int iColumn; for (iBlock = 0; iBlock < numberBlocks; iBlock++) { int nChunk; int jBlock; int iDo = iBlock * BLOCK; int base = iDo; if (iDo + BLOCK > numberRows_) { nChunk = numberRows_ - iDo; } else { nChunk = BLOCK; } solveF1(aa, nChunk, region + iDo); for (jBlock = iBlock + 1; jBlock < numberBlocks; jBlock++) { base += BLOCK; aa += BLOCKSQ; if (base + BLOCK > numberRows_) { nChunk = numberRows_ - base; } else { nChunk = BLOCK; } solveF2(aa, nChunk, region + iDo, region + base); } aa += BLOCKSQ; } /* do diagonal outside*/ for (iColumn = 0; iColumn < numberRows_; iColumn++) region[iColumn] *= diagonal_[iColumn]; int offset = ((numberBlocks * (numberBlocks + 1)) >> 1); aa = a + number_entries(offset - 1); int lBase = (numberBlocks - 1) * BLOCK; for (iBlock = numberBlocks - 1; iBlock >= 0; iBlock--) { int nChunk; int jBlock; int triBase = iBlock * BLOCK; int iBase = lBase; for (jBlock = iBlock + 1; jBlock < numberBlocks; jBlock++) { if (iBase + BLOCK > numberRows_) { nChunk = numberRows_ - iBase; } else { nChunk = BLOCK; } solveB2(aa, nChunk, region + triBase, region + iBase); iBase -= BLOCK; aa -= BLOCKSQ; } if (triBase + BLOCK > numberRows_) { nChunk = numberRows_ - triBase; } else { nChunk = BLOCK; } solveB1(aa, nChunk, region + triBase); aa -= BLOCKSQ; } #ifdef CHOL_COMPARE if (numberRows_ < 200) { for (int i = 0; i < numberRows_; i++) { assert(CoinAbs(region[i] - region2[i]) < 1.0e-3); } delete [] region2; } #endif } /* Forward part of solve 1*/ void ClpCholeskyDense::solveF1(longDouble * a, int n, CoinWorkDouble * region) { int j, k; CoinWorkDouble t00; for (j = 0; j < n; j ++) { t00 = region[j]; for (k = 0; k < j; ++k) { t00 -= region[k] * a[j + k * BLOCK]; } /*t00*=a[j + j * BLOCK];*/ region[j] = t00; } } /* Forward part of solve 2*/ void ClpCholeskyDense::solveF2(longDouble * a, int n, CoinWorkDouble * region, CoinWorkDouble * region2) { int j, k; #ifdef BLOCKUNROLL if (n == BLOCK) { for (k = 0; k < BLOCK; k += 4) { CoinWorkDouble t0 = region2[0]; CoinWorkDouble t1 = region2[1]; CoinWorkDouble t2 = region2[2]; CoinWorkDouble t3 = region2[3]; t0 -= region[0] * a[0 + 0 * BLOCK]; t1 -= region[0] * a[1 + 0 * BLOCK]; t2 -= region[0] * a[2 + 0 * BLOCK]; t3 -= region[0] * a[3 + 0 * BLOCK]; t0 -= region[1] * a[0 + 1 * BLOCK]; t1 -= region[1] * a[1 + 1 * BLOCK]; t2 -= region[1] * a[2 + 1 * BLOCK]; t3 -= region[1] * a[3 + 1 * BLOCK]; t0 -= region[2] * a[0 + 2 * BLOCK]; t1 -= region[2] * a[1 + 2 * BLOCK]; t2 -= region[2] * a[2 + 2 * BLOCK]; t3 -= region[2] * a[3 + 2 * BLOCK]; t0 -= region[3] * a[0 + 3 * BLOCK]; t1 -= region[3] * a[1 + 3 * BLOCK]; t2 -= region[3] * a[2 + 3 * BLOCK]; t3 -= region[3] * a[3 + 3 * BLOCK]; t0 -= region[4] * a[0 + 4 * BLOCK]; t1 -= region[4] * a[1 + 4 * BLOCK]; t2 -= region[4] * a[2 + 4 * BLOCK]; t3 -= region[4] * a[3 + 4 * BLOCK]; t0 -= region[5] * a[0 + 5 * BLOCK]; t1 -= region[5] * a[1 + 5 * BLOCK]; t2 -= region[5] * a[2 + 5 * BLOCK]; t3 -= region[5] * a[3 + 5 * BLOCK]; t0 -= region[6] * a[0 + 6 * BLOCK]; t1 -= region[6] * a[1 + 6 * BLOCK]; t2 -= region[6] * a[2 + 6 * BLOCK]; t3 -= region[6] * a[3 + 6 * BLOCK]; t0 -= region[7] * a[0 + 7 * BLOCK]; t1 -= region[7] * a[1 + 7 * BLOCK]; t2 -= region[7] * a[2 + 7 * BLOCK]; t3 -= region[7] * a[3 + 7 * BLOCK]; #if BLOCK>8 t0 -= region[8] * a[0 + 8 * BLOCK]; t1 -= region[8] * a[1 + 8 * BLOCK]; t2 -= region[8] * a[2 + 8 * BLOCK]; t3 -= region[8] * a[3 + 8 * BLOCK]; t0 -= region[9] * a[0 + 9 * BLOCK]; t1 -= region[9] * a[1 + 9 * BLOCK]; t2 -= region[9] * a[2 + 9 * BLOCK]; t3 -= region[9] * a[3 + 9 * BLOCK]; t0 -= region[10] * a[0 + 10 * BLOCK]; t1 -= region[10] * a[1 + 10 * BLOCK]; t2 -= region[10] * a[2 + 10 * BLOCK]; t3 -= region[10] * a[3 + 10 * BLOCK]; t0 -= region[11] * a[0 + 11 * BLOCK]; t1 -= region[11] * a[1 + 11 * BLOCK]; t2 -= region[11] * a[2 + 11 * BLOCK]; t3 -= region[11] * a[3 + 11 * BLOCK]; t0 -= region[12] * a[0 + 12 * BLOCK]; t1 -= region[12] * a[1 + 12 * BLOCK]; t2 -= region[12] * a[2 + 12 * BLOCK]; t3 -= region[12] * a[3 + 12 * BLOCK]; t0 -= region[13] * a[0 + 13 * BLOCK]; t1 -= region[13] * a[1 + 13 * BLOCK]; t2 -= region[13] * a[2 + 13 * BLOCK]; t3 -= region[13] * a[3 + 13 * BLOCK]; t0 -= region[14] * a[0 + 14 * BLOCK]; t1 -= region[14] * a[1 + 14 * BLOCK]; t2 -= region[14] * a[2 + 14 * BLOCK]; t3 -= region[14] * a[3 + 14 * BLOCK]; t0 -= region[15] * a[0 + 15 * BLOCK]; t1 -= region[15] * a[1 + 15 * BLOCK]; t2 -= region[15] * a[2 + 15 * BLOCK]; t3 -= region[15] * a[3 + 15 * BLOCK]; #endif region2[0] = t0; region2[1] = t1; region2[2] = t2; region2[3] = t3; region2 += 4; a += 4; } } else { #endif for (k = 0; k < n; ++k) { CoinWorkDouble t00 = region2[k]; for (j = 0; j < BLOCK; j ++) { t00 -= region[j] * a[k + j * BLOCK]; } region2[k] = t00; } #ifdef BLOCKUNROLL } #endif } /* Backward part of solve 1*/ void ClpCholeskyDense::solveB1(longDouble * a, int n, CoinWorkDouble * region) { int j, k; CoinWorkDouble t00; for (j = n - 1; j >= 0; j --) { t00 = region[j]; for (k = j + 1; k < n; ++k) { t00 -= region[k] * a[k + j * BLOCK]; } /*t00*=a[j + j * BLOCK];*/ region[j] = t00; } } /* Backward part of solve 2*/ void ClpCholeskyDense::solveB2(longDouble * a, int n, CoinWorkDouble * region, CoinWorkDouble * region2) { int j, k; #ifdef BLOCKUNROLL if (n == BLOCK) { for (j = 0; j < BLOCK; j += 4) { CoinWorkDouble t0 = region[0]; CoinWorkDouble t1 = region[1]; CoinWorkDouble t2 = region[2]; CoinWorkDouble t3 = region[3]; t0 -= region2[0] * a[0 + 0*BLOCK]; t1 -= region2[0] * a[0 + 1*BLOCK]; t2 -= region2[0] * a[0 + 2*BLOCK]; t3 -= region2[0] * a[0 + 3*BLOCK]; t0 -= region2[1] * a[1 + 0*BLOCK]; t1 -= region2[1] * a[1 + 1*BLOCK]; t2 -= region2[1] * a[1 + 2*BLOCK]; t3 -= region2[1] * a[1 + 3*BLOCK]; t0 -= region2[2] * a[2 + 0*BLOCK]; t1 -= region2[2] * a[2 + 1*BLOCK]; t2 -= region2[2] * a[2 + 2*BLOCK]; t3 -= region2[2] * a[2 + 3*BLOCK]; t0 -= region2[3] * a[3 + 0*BLOCK]; t1 -= region2[3] * a[3 + 1*BLOCK]; t2 -= region2[3] * a[3 + 2*BLOCK]; t3 -= region2[3] * a[3 + 3*BLOCK]; t0 -= region2[4] * a[4 + 0*BLOCK]; t1 -= region2[4] * a[4 + 1*BLOCK]; t2 -= region2[4] * a[4 + 2*BLOCK]; t3 -= region2[4] * a[4 + 3*BLOCK]; t0 -= region2[5] * a[5 + 0*BLOCK]; t1 -= region2[5] * a[5 + 1*BLOCK]; t2 -= region2[5] * a[5 + 2*BLOCK]; t3 -= region2[5] * a[5 + 3*BLOCK]; t0 -= region2[6] * a[6 + 0*BLOCK]; t1 -= region2[6] * a[6 + 1*BLOCK]; t2 -= region2[6] * a[6 + 2*BLOCK]; t3 -= region2[6] * a[6 + 3*BLOCK]; t0 -= region2[7] * a[7 + 0*BLOCK]; t1 -= region2[7] * a[7 + 1*BLOCK]; t2 -= region2[7] * a[7 + 2*BLOCK]; t3 -= region2[7] * a[7 + 3*BLOCK]; #if BLOCK>8 t0 -= region2[8] * a[8 + 0*BLOCK]; t1 -= region2[8] * a[8 + 1*BLOCK]; t2 -= region2[8] * a[8 + 2*BLOCK]; t3 -= region2[8] * a[8 + 3*BLOCK]; t0 -= region2[9] * a[9 + 0*BLOCK]; t1 -= region2[9] * a[9 + 1*BLOCK]; t2 -= region2[9] * a[9 + 2*BLOCK]; t3 -= region2[9] * a[9 + 3*BLOCK]; t0 -= region2[10] * a[10 + 0*BLOCK]; t1 -= region2[10] * a[10 + 1*BLOCK]; t2 -= region2[10] * a[10 + 2*BLOCK]; t3 -= region2[10] * a[10 + 3*BLOCK]; t0 -= region2[11] * a[11 + 0*BLOCK]; t1 -= region2[11] * a[11 + 1*BLOCK]; t2 -= region2[11] * a[11 + 2*BLOCK]; t3 -= region2[11] * a[11 + 3*BLOCK]; t0 -= region2[12] * a[12 + 0*BLOCK]; t1 -= region2[12] * a[12 + 1*BLOCK]; t2 -= region2[12] * a[12 + 2*BLOCK]; t3 -= region2[12] * a[12 + 3*BLOCK]; t0 -= region2[13] * a[13 + 0*BLOCK]; t1 -= region2[13] * a[13 + 1*BLOCK]; t2 -= region2[13] * a[13 + 2*BLOCK]; t3 -= region2[13] * a[13 + 3*BLOCK]; t0 -= region2[14] * a[14 + 0*BLOCK]; t1 -= region2[14] * a[14 + 1*BLOCK]; t2 -= region2[14] * a[14 + 2*BLOCK]; t3 -= region2[14] * a[14 + 3*BLOCK]; t0 -= region2[15] * a[15 + 0*BLOCK]; t1 -= region2[15] * a[15 + 1*BLOCK]; t2 -= region2[15] * a[15 + 2*BLOCK]; t3 -= region2[15] * a[15 + 3*BLOCK]; #endif region[0] = t0; region[1] = t1; region[2] = t2; region[3] = t3; a += 4 * BLOCK; region += 4; } } else { #endif for (j = 0; j < BLOCK; j ++) { CoinWorkDouble t00 = region[j]; for (k = 0; k < n; ++k) { t00 -= region2[k] * a[k + j * BLOCK]; } region[j] = t00; } #ifdef BLOCKUNROLL } #endif }