Changeset 423


Ignore:
Timestamp:
Sep 3, 2004 6:16:53 AM (15 years ago)
Author:
forrest
Message:

various Cholesky

Location:
trunk
Files:
12 edited

Legend:

Unmodified
Added
Removed
  • trunk/ClpCholeskyBase.cpp

    r414 r423  
    2525  type_(0),
    2626  doKKT_(false),
    27   zeroTolerance_(1.0e-17),
     27  goDense_(0.7),
    2828  choleskyCondition_(0.0),
    2929  model_(NULL),
     
    6363  type_(rhs.type_),
    6464  doKKT_(rhs.doKKT_),
    65   zeroTolerance_(rhs.zeroTolerance_),
     65  goDense_(rhs.goDense_),
    6666  choleskyCondition_(rhs.choleskyCondition_),
    6767  model_(rhs.model_),
     
    8282  choleskyRow_ = ClpCopyOfArray(rhs.choleskyRow_,sizeIndex_);
    8383  diagonal_ = ClpCopyOfArray(rhs.diagonal_,numberRows_);
     84#if CLP_LONG_CHOLESKY!=1
    8485  workDouble_ = ClpCopyOfArray(rhs.workDouble_,numberRows_);
     86#else
     87  // actually long double
     88  workDouble_ = (double *) ClpCopyOfArray((longWork *) rhs.workDouble_,numberRows_);
     89#endif
    8590  link_ = ClpCopyOfArray(rhs.link_,numberRows_);
    8691  workInteger_ = ClpCopyOfArray(rhs.workInteger_,numberRows_);
     
    127132    type_ = rhs.type_;
    128133    doKKT_ = rhs.doKKT_;
    129     zeroTolerance_ = rhs.zeroTolerance_;
     134    goDense_ = rhs.goDense_;
    130135    choleskyCondition_ = rhs.choleskyCondition_;
    131136    model_ = rhs.model_;
     
    162167    choleskyRow_ = ClpCopyOfArray(rhs.choleskyRow_,sizeIndex_);
    163168    diagonal_ = ClpCopyOfArray(rhs.diagonal_,numberRows_);
     169#if CLP_LONG_CHOLESKY!=1
    164170    workDouble_ = ClpCopyOfArray(rhs.workDouble_,numberRows_);
     171#else
     172    // actually long double
     173    workDouble_ = (double *) ClpCopyOfArray((longWork *) rhs.workDouble_,numberRows_);
     174#endif
    165175    link_ = ClpCopyOfArray(rhs.link_,numberRows_);
    166176    workInteger_ = ClpCopyOfArray(rhs.workInteger_,numberRows_);
     
    328338      } else {
    329339        // space for dense columns
    330         denseColumn_ = new double [numberDense*numberRows_];
     340        denseColumn_ = new longDouble [numberDense*numberRows_];
    331341        // dense cholesky
    332342        dense_ = new ClpCholeskyDense();
     
    623633  }
    624634  rowsDropped_ = new char [numberRows_];
    625   // use for counts
    626   memset(rowsDropped_,0,numberRows_);
    627635  numberRowsDropped_=0;
    628636  rowCopy_ = model->clpMatrix()->reverseOrderedCopy();
     
    633641  const int * rowLength = rowCopy_->getVectorLengths();
    634642  const int * column = rowCopy_->getIndices();
    635   // We need two arrays for counts
     643  // We need arrays for counts
    636644  int * which = new int [numberRows_];
    637645  int * used = new int[numberRows_+1];
     646  int * count = new int[numberRows_];
     647  CoinZeroN(count,numberRows_);
    638648  CoinZeroN(used,numberRows_);
    639649  int iRow;
     
    681691      } else {
    682692        // space for dense columns
    683         denseColumn_ = new double [numberDense*numberRows_];
     693        denseColumn_ = new longDouble [numberDense*numberRows_];
    684694        // dense cholesky
    685695        dense_ = new ClpCholeskyDense();
     
    709719                used[jRow]=1;
    710720                which[number++]=jRow;
    711                 rowsDropped_[jRow]++;
     721                count[jRow]++;
    712722              }
    713723            }
    714724          }
    715725        }
    716     }
     726      }
    717727      sizeFactor_ += number;
    718       rowsDropped_[iRow]+=number;
     728      count[iRow]+=number;
    719729      int j;
    720730      for (j=0;j<number;j++)
    721731        used[which[j]]=0;
    722732    }
    723     CoinSort_2(rowsDropped_,rowsDropped_+numberRows_,permute_);
     733    CoinSort_2(count,count+numberRows_,permute_);
    724734  } else {
    725735    // KKT
     
    735745  delete [] which;
    736746  delete [] used;
     747  delete [] count;
    737748  permuteInverse_ = new int [numberRows_];
    738749  memset(rowsDropped_,0,numberRows_);
     
    745756  return 0;
    746757}
    747 static int yyyyyy=0;
    748758/* Does Symbolic factorization given permutation.
    749759   This is called immediately after order.  If user provides this then
     
    769779  // We need an array for counts
    770780  int * used = new int[numberRows_+1];
     781  // If KKT then re-order so negative first
     782  if (doKKT_) {
     783    int nn=0;
     784    int np=0;
     785    int iRow;
     786    for (iRow=0;iRow<numberRows_;iRow++) {
     787      int originalRow = permute_[iRow];
     788      if (originalRow<numberTotal)
     789        permute_[nn++]=originalRow;
     790      else
     791        used[np++]=originalRow;
     792    }
     793    memcpy(permute_+nn,used,np*sizeof(int));
     794    for (iRow=0;iRow<numberRows_;iRow++)
     795      permuteInverse_[permute_[iRow]]=iRow;
     796  }
    771797  CoinZeroN(used,numberRows_);
    772798  int iRow;
     
    839865      }
    840866    }
    841     // but fake it
    842     if (yyyyyy==1) {
    843       for (iRow=0;iRow<numberRows_;iRow++) {
    844         //permute_[iRow]=iRow; // force no permute
    845         permute_[iRow]=numberRows_-1-iRow; // force odd permute
    846         //permute_[iRow]=(iRow+1)%numberRows_; // force odd permute
    847         permuteInverse_[permute_[iRow]]=iRow;
    848       }
    849       permuted=true;
    850     } else if (yyyyyy==2) {
    851       for (iRow=0;iRow<numberRows_;iRow++) {
    852         permute_[iRow]=iRow; // force no permute
    853         //permute_[iRow]=numberRows_-1-iRow; // force odd permute
    854         //permute_[iRow]=(iRow+1)%numberRows_; // force odd permute
    855         permuteInverse_[permute_[iRow]]=iRow;
    856       }
    857     }
    858867    if (permuted) {
    859868      // Need to permute - ugly
     
    10041013    noMemory=true;
    10051014  }
     1015  double sizeFactor=sizeFactor_;
    10061016  if (!noMemory) {
    10071017    // Do lower triangular
     
    12161226  }
    12171227  if (model_->messageHandler()->logLevel()>0)
    1218     std::cout<<sizeFactor_<<" elements in sparse Cholesky, flop count "<<flops<<std::endl;
     1228    std::cout<<sizeFactor<<" elements in sparse Cholesky, flop count "<<flops<<std::endl;
    12191229  try {
    1220     sparseFactor_ = new double [sizeFactor_];
    1221     workDouble_ = new double[numberRows_];
    1222     diagonal_ = new double[numberRows_];
     1230    sparseFactor_ = new longDouble [sizeFactor_];
     1231#if CLP_LONG_CHOLESKY!=1
     1232    workDouble_ = new longDouble[numberRows_];
     1233#else
     1234    // actually long double
     1235    workDouble_ = (double *) (new longWork[numberRows_]);
     1236#endif
     1237    diagonal_ = new longDouble[numberRows_];
    12231238  }
    12241239  catch (...) {
     
    12551270  return 0;
    12561271}
    1257 static int xxxxxx=2;
    12581272int
    12591273ClpCholeskyBase::symbolic1(const CoinBigIndex * Astart, const int * Arow)
    12601274{
    1261   if (xxxxxx!=2) {
    1262     int iRow;
    1263     int put=0;
    1264     choleskyStart_[0]=put;
    1265     for (iRow=0;iRow<numberRows_;iRow++) {
    1266       put += numberRows_-iRow-1;
    1267       choleskyStart_[iRow+1]=put;
    1268       clique_[iRow]=-1;
    1269     }
    1270     sizeFactor_ = choleskyStart_[numberRows_];
    1271     return sizeFactor_;
    1272   }
    12731275  int * marked = (int *) workInteger_;
    12741276  int iRow;
     
    13041306ClpCholeskyBase::symbolic2(const CoinBigIndex * Astart, const int * Arow)
    13051307{
    1306   if (xxxxxx!=2) {
    1307     int iRow;
    1308     int put=0;
    1309     choleskyStart_[0]=put;
    1310     for (iRow=0;iRow<numberRows_;iRow++) {
    1311       put += numberRows_-iRow-1;
    1312       choleskyStart_[iRow+1]=put;
    1313       clique_[iRow]=-1;
    1314     }
    1315     memcpy(indexStart_,choleskyStart_,numberRows_*sizeof(int));
    1316     sizeFactor_ = choleskyStart_[numberRows_];
    1317     sizeIndex_ = sizeFactor_;
    1318     put=0;
    1319    
    1320     xxxxxx=2;
    1321     for (iRow=0;iRow<numberRows_;iRow++) {
    1322       for (int j=iRow+1;j<numberRows_;j++) {
    1323         choleskyRow_[put++]=j;
    1324       }
    1325     }
    1326     return;
    1327   }
    1328   int iRow;
    1329 #if 0
    1330   int * marked = (int *) workInteger_;
    1331   int * first = new int[numberRows_];
    1332   int * index = new int[numberRows_];
    1333   sizeFactor_=0;
    1334   sizeIndex_=0;
    1335   choleskyStart_[0]=0;
    1336   // may not need to do this here but makes debugging easier
    1337   for (iRow=0;iRow<numberRows_;iRow++) {
    1338     marked[iRow]=-1;
    1339     link_[iRow]=-1;
    1340     first[iRow]=Astart[iRow];
    1341   }
    1342   for (iRow=0;iRow<numberRows_;iRow++) {
    1343     int kRow=iRow;
    1344     int number=0;
    1345     // skip diagonal
    1346     int n=Astart[iRow+1]-Astart[iRow]-1;
    1347     choleskyStart_[iRow+1]=sizeFactor_+n;
    1348     memcpy(choleskyRow_+sizeFactor_,Arow+Astart[iRow]+1,n*sizeof(int));
    1349     first[iRow]=sizeFactor_;
    1350     while (kRow>=0) {
    1351       int nextRow=link_[kRow];
    1352       CoinBigIndex k=first[kRow];
    1353       CoinBigIndex end=choleskyStart_[kRow+1];
    1354       if (k+1<end) {
    1355         // update first
    1356         if (kRow!=iRow) {
    1357           assert (iRow==choleskyRow_[k]);
    1358           k++;
    1359           first[kRow]=k;
    1360         }
    1361         int jRow = choleskyRow_[k];
    1362         link_[kRow]=link_[jRow];
    1363         link_[jRow]=kRow;
    1364         for (;k<end;k++) {
    1365           int jRow = choleskyRow_[k];
    1366           if (marked[jRow]<0) {
    1367             // not marked yet
    1368             // may want to compare indices to last (later)
    1369             marked[jRow]=1;
    1370             index[number++]=jRow;
    1371           }
    1372         }
    1373       }
    1374       kRow=nextRow;
    1375     }
    1376     CoinBigIndex put=sizeFactor_;
    1377     // could look for clique
    1378     for (int i=0;i<number;i++) {
    1379       int jRow=index[i];
    1380       assert (jRow!=iRow);
    1381       choleskyRow_[put++]=jRow;
    1382       marked[jRow]=-1;
    1383     }
    1384     sizeFactor_ = put;
    1385     choleskyStart_[iRow+1]=sizeFactor_;
    1386     indexStart_[iRow]=sizeIndex_;
    1387     sizeIndex_ =put;
    1388   }
    1389   delete [] first;
    1390   delete [] index;
    1391   for (iRow=0;iRow<numberRows_;iRow++)
    1392     clique_[iRow]=-1;
    1393   return;
    1394 #endif
    13951308  int * mergeLink = clique_;
    13961309  int * marker = (int *) workInteger_;
    1397   //int iRow;
     1310  int iRow;
    13981311  for (iRow=0;iRow<numberRows_;iRow++) {
    13991312    marker[iRow]=-1;
     
    14321345    bool reuse=false;
    14331346    // Check if we can re-use indices
    1434     if (!marked && merge>=0&&mergeLink[merge]<0&&0) {
     1347    if (!marked && merge>=0&&mergeLink[merge]<0) {
    14351348      // can re-use all
    1436       indexStart_[iRow]=indexStart_[merge]+1;
     1349      startSub = indexStart_[merge]+1;
    14371350      nz=choleskyStart_[merge+1]-(choleskyStart_[merge]+1);
    14381351      reuse=true;
     
    14701383        reuse=true; // can re-use
    14711384    }
    1472     reuse=false; //temp
     1385    //reuse=false; //temp
    14731386    if (!reuse) {
    14741387      end += nz;
     
    14911404    }
    14921405    // should not be needed
    1493     std::sort(choleskyRow_+choleskyStart_[iRow]
    1494               ,choleskyRow_+choleskyStart_[iRow+1]);
     1406    //std::sort(choleskyRow_+indexStart_[iRow]
     1407    //      ,choleskyRow_+indexStart_[iRow]+nz);
     1408#define CLP_DEBUG
    14951409#ifdef CLP_DEBUG
    14961410    int last=-1;
    1497     for (CoinBigIndex j=indexStart_[k];j<startSub;j++) {
     1411    for (CoinBigIndex j=indexStart_[iRow];j<indexStart_[iRow]+nz;j++) {
    14981412      int kRow=choleskyRow_[j];
    14991413      assert (kRow>last);
     
    15051419  sizeIndex_ = start;
    15061420  // find dense segment here
     1421  int numberleft=numberRows_;
     1422  for (iRow=0;iRow<numberRows_;iRow++) {
     1423    CoinBigIndex left=sizeFactor_-choleskyStart_[iRow];
     1424    double n=numberleft;
     1425    double threshold = n*(n-1.0)*0.5*goDense_;
     1426    if ((double) left >= threshold)
     1427      break;
     1428    numberleft--;
     1429  }
     1430  //iRow=numberRows_;
     1431  int nDense = numberRows_-iRow;
     1432#define DENSE_THRESHOLD 8
     1433  // don't do if dense columns
     1434  if (nDense>=DENSE_THRESHOLD&&!dense_) {
     1435    printf("Going dense for last %d rows\n",nDense);
     1436    // make sure we don't disturb any indices
     1437    CoinBigIndex k=0;
     1438    for (int jRow=0;jRow<iRow;jRow++) {
     1439      int nz=choleskyStart_[jRow+1]-choleskyStart_[jRow];
     1440      k=CoinMax(k,indexStart_[jRow]+nz);
     1441    }
     1442    indexStart_[iRow]=k;
     1443    int j;
     1444    for (j=iRow+1;j<numberRows_;j++) {
     1445      choleskyRow_[k++]=j;
     1446      indexStart_[j]=k;
     1447    }
     1448    sizeIndex_=k;
     1449    assert (k<=sizeFactor_); // can't happen with any reasonable defaults
     1450    k=choleskyStart_[iRow];
     1451    for (j=iRow+1;j<=numberRows_;j++) {
     1452      k += numberRows_-j;
     1453      choleskyStart_[j]=k;
     1454    }
     1455    // allow for blocked dense
     1456    ClpCholeskyDense dense;
     1457    sizeFactor_=choleskyStart_[iRow]+dense.space(nDense);
     1458    firstDense_=iRow;
     1459    if (doKKT_) {
     1460      // redo permute so negative ones first
     1461      int putN=firstDense_;
     1462      int putP=0;
     1463      int numberRowsModel = model_->numberRows();
     1464      int numberColumns = model_->numberColumns();
     1465      int numberTotal = numberColumns + numberRowsModel;
     1466      for (iRow=firstDense_;iRow<numberRows_;iRow++) {
     1467        int originalRow=permute_[iRow];
     1468        if (originalRow<numberTotal)
     1469          permute_[putN++]=originalRow;
     1470        else
     1471          permuteInverse_[putP++]=originalRow;
     1472      }
     1473      for (iRow=putN;iRow<numberRows_;iRow++) {
     1474        permute_[iRow]=permuteInverse_[iRow-putN];
     1475      }
     1476      for (iRow=0;iRow<numberRows_;iRow++) {
     1477        permuteInverse_[permute_[iRow]]=iRow;
     1478      }
     1479    }
     1480  }
    15071481  // Clean up clique info
    15081482  for (iRow=0;iRow<numberRows_;iRow++)
    1509     clique_[iRow]=-1;
    1510   // should not be needed
    1511   for (iRow=0;iRow<numberRows_;iRow++) {
    1512     std::sort(choleskyRow_+choleskyStart_[iRow],
    1513               choleskyRow_+choleskyStart_[iRow+1]);
     1483    clique_[iRow]=0;
     1484  int lastClique=-1;
     1485  bool inClique=false;
     1486  for (iRow=1;iRow<firstDense_;iRow++) {
     1487    int sizeLast = choleskyStart_[iRow]-choleskyStart_[iRow-1];
     1488    int sizeThis = choleskyStart_[iRow+1]-choleskyStart_[iRow];
     1489    if (indexStart_[iRow]==indexStart_[iRow-1]+1&&
     1490        sizeThis==sizeLast-1) {
     1491      // in clique
     1492      if (!inClique) {
     1493        inClique=true;
     1494        lastClique=iRow-1;
     1495      }
     1496    } else if (inClique) {
     1497      int sizeClique=iRow-lastClique;
     1498      for (int i=lastClique;i<iRow;i++) {
     1499        clique_[i]=sizeClique;
     1500        sizeClique--;
     1501      }
     1502      inClique=false;
     1503    }
     1504  }
     1505  if (inClique) {
     1506    int sizeClique=iRow-lastClique;
     1507    for (int i=lastClique;i<iRow;i++) {
     1508      clique_[i]=sizeClique;
     1509      sizeClique--;
     1510    }
    15141511  }
    15151512}
     
    15281525  int numberColumns=model_->clpMatrix()->getNumCols();
    15291526  //perturbation
    1530   double perturbation=model_->diagonalPerturbation()*model_->diagonalNorm();
     1527  longDouble perturbation=model_->diagonalPerturbation()*model_->diagonalNorm();
    15311528  perturbation=perturbation*perturbation;
    15321529  if (perturbation>1.0) {
     
    15381535  int iRow;
    15391536  int iColumn;
    1540   double * work = workDouble_;
     1537  longDouble * work = workDouble_;
    15411538  CoinZeroN(work,numberRows_);
    15421539  int newDropped=0;
     
    15491546      numberDense=dense_->numberRows();
    15501547    if (whichDense_) {
    1551       double * denseDiagonal = dense_->diagonal();
    1552       double * dense = denseColumn_;
     1548      longDouble * denseDiagonal = dense_->diagonal();
     1549      longDouble * dense = denseColumn_;
    15531550      int iDense=0;
    15541551      for (int iColumn=0;iColumn<numberColumns;iColumn++) {
    15551552        if (whichDense_[iColumn]) {
    1556           denseDiagonal[iDense++]=1.0/diagonal[iColumn];
    15571553          CoinZeroN(dense,numberRows_);
    15581554          CoinBigIndex start=columnStart[iColumn];
    15591555          CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
    1560           for (CoinBigIndex j=start;j<end;j++) {
    1561             int jRow=row[j];
    1562             dense[jRow] = element[j];
     1556          if (diagonal[iColumn]) {
     1557            denseDiagonal[iDense++]=1.0/diagonal[iColumn];
     1558            for (CoinBigIndex j=start;j<end;j++) {
     1559              int jRow=row[j];
     1560              dense[jRow] = element[j];
     1561            }
     1562          } else {
     1563            denseDiagonal[iDense++]=1.0;
    15631564          }
    15641565          dense += numberRows_;
     
    15661567      }
    15671568    }
    1568     double delta2 = model_->delta(); // add delta*delta to diagonal
     1569    longDouble delta2 = model_->delta(); // add delta*delta to diagonal
    15691570    delta2 *= delta2;
    15701571    // largest in initial matrix
    15711572    double largest2=1.0e-20;
    15721573    for (iRow=0;iRow<numberRows_;iRow++) {
    1573       double * put = sparseFactor_+choleskyStart_[iRow];
     1574      longDouble * put = sparseFactor_+choleskyStart_[iRow];
    15741575      int * which = choleskyRow_+indexStart_[iRow];
    15751576      int iOriginalRow = permute_[iRow];
     
    15861587            CoinBigIndex start=columnStart[iColumn];
    15871588            CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
    1588             double multiplier = diagonal[iColumn]*elementByRow[k];
     1589            longDouble multiplier = diagonal[iColumn]*elementByRow[k];
    15891590            for (CoinBigIndex j=start;j<end;j++) {
    15901591              int jRow=row[j];
    15911592              int jNewRow = permuteInverse_[jRow];
    15921593              if (jNewRow>=iRow&&!rowsDropped_[jRow]) {
    1593                 double value=element[j]*multiplier;
     1594                longDouble value=element[j]*multiplier;
    15941595                work[jNewRow] += value;
    15951596              }
     
    16251626      rowsDropped[iRow]=dropped;
    16261627      if (!dropped) {
    1627         double diagonal = diagonal_[iRow];
     1628        longDouble diagonal = diagonal_[iRow];
    16281629        if (diagonal>largest2) {
    16291630          diagonal_[iRow]=diagonal+perturbation;
     
    16561657      int i;
    16571658      for (i=0;i<numberDense;i++) {
    1658         double * a = denseColumn_+i*numberRows_;
     1659        longDouble * a = denseColumn_+i*numberRows_;
    16591660        for (int j=0;j<numberRows_;j++) {
    16601661          if (rowsDropped[j])
    16611662            a[j]=0.0;
    16621663        }
    1663         solve(a,1);
     1664        solveLong(a,1);
    16641665      }
    16651666      dense_->resetRowsDropped();
    16661667      longDouble * denseBlob = dense_->aMatrix();
    1667       double * denseDiagonal = dense_->diagonal();
     1668      longDouble * denseDiagonal = dense_->diagonal();
    16681669      // Update dense matrix
    16691670      for (i=0;i<numberDense;i++) {
    1670         const double * a = denseColumn_+i*numberRows_;
     1671        const longDouble * a = denseColumn_+i*numberRows_;
    16711672        // do diagonal
    1672         double value = denseDiagonal[i];
    1673         const double * b = denseColumn_+i*numberRows_;
     1673        longDouble value = denseDiagonal[i];
     1674        const longDouble * b = denseColumn_+i*numberRows_;
    16741675        for (int k=0;k<numberRows_;k++)
    16751676          value += a[k]*b[k];
    16761677        denseDiagonal[i]=value;
    16771678        for (int j=i+1;j<numberDense;j++) {
    1678           double value = 0.0;
    1679           const double * b = denseColumn_+j*numberRows_;
     1679          longDouble value = 0.0;
     1680          const longDouble * b = denseColumn_+j*numberRows_;
    16801681          for (int k=0;k<numberRows_;k++)
    16811682            value += a[k]*b[k];
     
    16891690      delete [] dropped;
    16901691    }
     1692    // try allowing all every time
     1693    //printf("trying ?\n");
     1694    //for (iRow=0;iRow<numberRows_;iRow++) {
     1695    //rowsDropped[iRow]=0;
     1696    //rowsDropped_[iRow]=0;
     1697    //}
    16911698    bool cleanCholesky;
     1699    //if (model_->numberIterations()<20||(model_->numberIterations()&1)==0)
    16921700    if (model_->numberIterations()<2000)
    16931701      cleanCholesky=true;
     
    17011709          char dropped = rowsDropped[i];
    17021710          rowsDropped_[i]=dropped;
     1711          rowsDropped_[i]=0;
    17031712          if (dropped==2) {
    17041713            //dropped this time
     
    17581767    }
    17591768    if (permuted) {
    1760       double delta2 = model_->delta(); // add delta*delta to bottom
     1769      longDouble delta2 = model_->delta(); // add delta*delta to bottom
    17611770      delta2 *= delta2;
    17621771      // Need to permute - ugly
    17631772      if (!quadratic) {
    17641773        for (iRow=0;iRow<numberRows_;iRow++) {
    1765           double * put = sparseFactor_+choleskyStart_[iRow];
     1774          longDouble * put = sparseFactor_+choleskyStart_[iRow];
    17661775          int * which = choleskyRow_+indexStart_[iRow];
    17671776          int iOriginalRow = permute_[iRow];
    17681777          if (iOriginalRow<numberColumns) {
    17691778            iColumn=iOriginalRow;
    1770             double value = diagonal[iColumn];
     1779            longDouble value = diagonal[iColumn];
    17711780            if (fabs(value)>1.0e-100) {
    17721781              value = 1.0/value;
     
    17871796            }
    17881797          } else if (iOriginalRow<numberTotal) {
    1789             double value = diagonal[iOriginalRow];
     1798            longDouble value = diagonal[iOriginalRow];
    17901799            if (fabs(value)>1.0e-100) {
    17911800              value = 1.0/value;
     
    18311840        const double * quadraticElement = quadratic->getElements();
    18321841        for (iRow=0;iRow<numberRows_;iRow++) {
    1833           double * put = sparseFactor_+choleskyStart_[iRow];
     1842          longDouble * put = sparseFactor_+choleskyStart_[iRow];
    18341843          int * which = choleskyRow_+indexStart_[iRow];
    18351844          int iOriginalRow = permute_[iRow];
     
    18371846            CoinBigIndex j;
    18381847            iColumn=iOriginalRow;
    1839             double value = diagonal[iColumn];
     1848            longDouble value = diagonal[iColumn];
    18401849            if (fabs(value)>1.0e-100) {
    18411850              value = 1.0/value;
     
    18661875            }
    18671876          } else if (iOriginalRow<numberTotal) {
    1868             double value = diagonal[iOriginalRow];
     1877            longDouble value = diagonal[iOriginalRow];
    18691878            if (fabs(value)>1.0e-100) {
    18701879              value = 1.0/value;
     
    19091918      if (!quadratic) {
    19101919        for (iColumn=0;iColumn<numberColumns;iColumn++) {
    1911           double * put = sparseFactor_+choleskyStart_[iColumn];
     1920          longDouble * put = sparseFactor_+choleskyStart_[iColumn];
    19121921          int * which = choleskyRow_+indexStart_[iColumn];
    1913           double value = diagonal[iColumn];
     1922          longDouble value = diagonal[iColumn];
    19141923          if (fabs(value)>1.0e-100) {
    19151924            value = 1.0/value;
     
    19421951        const double * quadraticElement = quadratic->getElements();
    19431952        for (iColumn=0;iColumn<numberColumns;iColumn++) {
    1944           double * put = sparseFactor_+choleskyStart_[iColumn];
     1953          longDouble * put = sparseFactor_+choleskyStart_[iColumn];
    19451954          int * which = choleskyRow_+indexStart_[iColumn];
    19461955          int number = choleskyStart_[iColumn+1]-choleskyStart_[iColumn];
    1947           double value = diagonal[iColumn];
     1956          longDouble value = diagonal[iColumn];
    19481957          CoinBigIndex j;
    19491958          if (fabs(value)>1.0e-100) {
     
    19831992      // slacks
    19841993      for (iColumn=numberColumns;iColumn<numberTotal;iColumn++) {
    1985         double * put = sparseFactor_+choleskyStart_[iColumn];
     1994        longDouble * put = sparseFactor_+choleskyStart_[iColumn];
    19861995        int * which = choleskyRow_+indexStart_[iColumn];
    1987         double value = diagonal[iColumn];
     1996        longDouble value = diagonal[iColumn];
    19881997        if (fabs(value)>1.0e-100) {
    19891998          value = 1.0/value;
     
    20032012      }
    20042013      // Finish diagonal
    2005       double delta2 = model_->delta(); // add delta*delta to bottom
     2014      longDouble delta2 = model_->delta(); // add delta*delta to bottom
    20062015      delta2 *= delta2;
    20072016      for (iRow=0;iRow<numberRowsModel;iRow++) {
    2008         double * put = sparseFactor_+choleskyStart_[iRow+numberTotal];
     2017        longDouble * put = sparseFactor_+choleskyStart_[iRow+numberTotal];
    20092018        diagonal_[iRow+numberTotal]=delta2;
    20102019        CoinBigIndex j;
     
    20662075  return newDropped;
    20672076}
    2068 static /* Subroutine */ int ekkagmyldlt3(double *a, int *lda, int *n,
    2069                                   double *v, double *din)
    2070 {
    2071     /* System generated locals */
    2072     int a_dim1, i__1, i__2;
    2073 
    2074     /* Local variables */
    2075     int i, j, k;
    2076     double y, x0;
    2077     double t00;
    2078     double ai0, aj0, rd0;
    2079     double yj0;
    2080     double ddmin, ddmax;
    2081     double dsmall;
    2082 
    2083 /*      A contains the original matrix and the result is returned in lvals. */
    2084 /*      A is garbage on output.  Both A and lvals are lower triangular normal.
    2085  */
    2086 /*     do j = 0, N - 1 */
    2087 /*       do i = 0, j - 1 */
    2088 /*         x = v(i) * A(j,i) */
    2089 /*         do k = j, N - 1 */
    2090 /*           A(k,j) = A(k,j) - A(k,i) * x */
    2091 /*         end do */
    2092 /*       end do */
    2093 /*       x = A(j,j) */
    2094 /*       v(j) = x */
    2095 /*       y = one/x */
    2096 /*       A(j,j) = y */
    2097 /*       do i = j + 1, N - 1 */
    2098 /*         A(i,j) = A(i,j) * y */
    2099 /*       end do */
    2100 /*       print *,'Matrix after scalar iter = ',j */
    2101 /*      call prtmat(A,lda,N) */
    2102 /*     end do */
    2103 /*     return */
    2104     /* Parameter adjustments */
    2105     a_dim1 = *lda;
    2106 
    2107     /* Function Body */
    2108     dsmall = 1.0e-20;
    2109     ddmax = 0.0;
    2110     ddmin = 1.0e30;
    2111     i__1 = *n;
    2112     for (j = 0; j < i__1; j ++) {
    2113 
    2114 /* process column j */
    2115 
    2116       //t00 = a[j + j * a_dim1];
    2117       t00 = din[j];
    2118       assert (t00==din[j]);
    2119       //assert (t00==v[j]);
    2120       for (k = 0; k <j; ++k) {
    2121         y = v[k];
    2122         aj0 = a[j + k * a_dim1];
    2123         yj0 = aj0 * y;
    2124         t00 -= aj0 * yj0;
    2125       }
    2126       x0 = fabs(t00);
    2127       if (x0 > ddmax) {
    2128         ddmax = t00;
    2129       }
    2130       if (x0 < ddmin) {
    2131         ddmin = t00;
    2132       }
    2133       if (x0 <= dsmall) {
    2134         abort();
    2135       }
    2136       rd0 = 1. / t00;
    2137       v[j] = t00;
    2138       //a[j + j * a_dim1] = rd0;
    2139       din[j]=rd0;
    2140 
    2141 /* process the remainder of column j */
    2142 
    2143       i__2 = *n;
    2144       for (i = j+1; i < i__2; i ++) {
    2145         t00 = a[i + j * a_dim1];
    2146         for (k = 0; k <j; ++k) {
    2147           y = v[k];
    2148           aj0 = a[j + k * a_dim1] * y;
    2149           ai0 = a[i + k * a_dim1];
    2150           t00 -= ai0 * aj0;
    2151         }
    2152         t00 *= rd0;
    2153         a[i + j * a_dim1] = t00;
    2154       }
    2155     }
    2156     return 0;
    2157 
    2158 } /* ekkagmyldlt3 */
    21592077/* Factorize - filling in rowsDropped and returning number dropped
    21602078   in integerParam.
     
    21632081ClpCholeskyBase::factorizePart2(int * rowsDropped)
    21642082{
    2165   double & largest=doubleParameters_[3];
    2166   double & smallest=doubleParameters_[4];
     2083  longDouble largest=doubleParameters_[3];
     2084  longDouble smallest=doubleParameters_[4];
    21672085  int numberDropped=integerParameters_[20];
    21682086  // probably done before
     
    21722090  double dropValue = doubleParameters_[10];
    21732091  int firstPositive=integerParameters_[34];
    2174   double * d = ClpCopyOfArray(diagonal_,numberRows_);
     2092  longDouble * d = ClpCopyOfArray(diagonal_,numberRows_);
    21752093  int iRow;
    2176   if (0) {
    2177     double aa[40000];
    2178     //double v[200];
    2179     assert (numberRows_<200);
    2180     memset(d,0,numberRows_*sizeof(double));
    2181     memset(aa,0,numberRows_*numberRows_*sizeof(double));
    2182     for (iRow=0;iRow<numberRows_;iRow++) {
    2183       int base = numberRows_*iRow;
    2184       aa[base+iRow]=diagonal_[iRow];
    2185       for (int k=choleskyStart_[iRow];k<choleskyStart_[iRow+1];k++) {
    2186         int jRow = choleskyRow_[k];
    2187         aa[base+jRow]=sparseFactor_[k];
    2188       }
    2189     }
    2190     ekkagmyldlt3(aa, &numberRows_,&numberRows_,d,diagonal_);
    2191     //memcpy(diagonal_,v,numberRows_*sizeof(double));
    2192     for (iRow=0;iRow<numberRows_;iRow++) {
    2193       int base = numberRows_*iRow;
    2194       for (int k=choleskyStart_[iRow];k<choleskyStart_[iRow+1];k++) {
    2195         int jRow = choleskyRow_[k];
    2196         sparseFactor_[k]=aa[base+jRow];
    2197       }
    2198     }
    2199     return;
    2200   }
    22012094  // minimum size before clique done
    22022095#define MINCLIQUE INT_MAX
    2203   double * work = workDouble_;
     2096  longDouble * work = workDouble_;
    22042097  CoinBigIndex * first = workInteger_;
    22052098 
     
    22142107  bool newClique=false;
    22152108  bool endClique=false;
     2109  int lastRow=0;
     2110  int cliqueSize=0;
     2111  CoinBigIndex cliquePointer=0;
     2112  int nextRow2=-1;
    22162113 
    2217   for (iRow=0;iRow<firstDense_;iRow++) {
    2218     endClique=false;
    2219     if (clique_[iRow]>=0) {
    2220       // this is in a clique
    2221       inClique=true;
    2222       if (clique_[iRow]>lastClique) {
    2223         // new Clique
    2224         newClique=true;
     2114  for (iRow=0;iRow<firstDense_+1;iRow++) {
     2115    if (iRow<firstDense_) {
     2116      endClique=false;
     2117      if (clique_[iRow]>0) {
     2118        // this is in a clique
     2119        inClique=true;
     2120        if (clique_[iRow]>lastClique) {
     2121          // new Clique
     2122          newClique=true;
     2123          // If we have clique going then signal to do old one
     2124          endClique=(lastClique>0);
     2125        } else {
     2126          // Still in clique
     2127          newClique=false;
     2128        }
     2129      } else {
     2130        // not in clique
     2131        inClique=false;
     2132        newClique=false;
    22252133        // If we have clique going then signal to do old one
    2226         endClique=(lastClique>=0);
    2227       } else {
    2228         // Still in clique
    2229         newClique=false;
    2230       }
     2134        endClique=(lastClique>0);
     2135      }
     2136      lastClique=clique_[iRow];
     2137    } else if (inClique) {
     2138      // Finish off
     2139      endClique=true;
    22312140    } else {
    2232       // not in clique
    2233       inClique=false;
    2234       newClique=false;
    2235       // If we have clique going then signal to do old one
    2236       endClique=(lastClique>=0);
    2237     }
    2238     lastClique=clique_[iRow];
     2141      break;
     2142    }
    22392143    if (endClique) {
    22402144      // We have just finished updating a clique - do block pivot and clean up
    2241       abort(); //ekkchpp coding
    2242     }
     2145      int jRow;
     2146      for ( jRow=lastRow;jRow<iRow;jRow++) {
     2147        int jCount = jRow-lastRow;
     2148        longDouble diagonalValue = diagonal_[jRow];
     2149        CoinBigIndex start=choleskyStart_[jRow];
     2150        CoinBigIndex end=choleskyStart_[jRow+1];
     2151        for (int kRow=lastRow;kRow<jRow;kRow++) {
     2152          jCount--;
     2153          CoinBigIndex get = choleskyStart_[kRow]+jCount;
     2154          longDouble a_jk = sparseFactor_[get];
     2155          longDouble value1 = d[kRow]*a_jk;
     2156          diagonalValue -= a_jk*value1;
     2157          for (CoinBigIndex j=start;j<end;j++)
     2158            sparseFactor_[j] -= value1*sparseFactor_[++get];
     2159        }
     2160        // check
     2161        int originalRow = permute_[jRow];
     2162        if (originalRow<firstPositive) {
     2163          // must be negative
     2164          if (diagonalValue<=-dropValue) {
     2165            smallest = CoinMin(smallest,-diagonalValue);
     2166            largest = CoinMax(largest,-diagonalValue);
     2167            d[jRow]=diagonalValue;
     2168            diagonalValue = 1.0/diagonalValue;
     2169          } else {
     2170            rowsDropped[originalRow]=2;
     2171            d[jRow]=-1.0e100;
     2172            diagonalValue=0.0;
     2173            integerParameters_[20]++;
     2174          }
     2175        } else {
     2176          // must be positive
     2177          if (diagonalValue>=dropValue) {
     2178            smallest = CoinMin(smallest,diagonalValue);
     2179            largest = CoinMax(largest,diagonalValue);
     2180            d[jRow]=diagonalValue;
     2181            diagonalValue = 1.0/diagonalValue;
     2182          } else {
     2183            rowsDropped[originalRow]=2;
     2184            d[jRow]=1.0e100;
     2185            diagonalValue=0.0;
     2186            integerParameters_[20]++;
     2187          }
     2188        }
     2189        diagonal_[jRow]=diagonalValue;
     2190        for (CoinBigIndex j=start;j<end;j++) {
     2191          sparseFactor_[j] *= diagonalValue;
     2192        }
     2193      }
     2194      if (nextRow2>=0) {
     2195        for ( jRow=lastRow;jRow<iRow-1;jRow++) {
     2196          link_[jRow]=jRow+1;
     2197        }
     2198        link_[iRow-1]=link_[nextRow2];
     2199        link_[nextRow2]=lastRow;
     2200      }
     2201    }
     2202    if (iRow==firstDense_)
     2203      break; // we were just cleaning up
    22432204    if (newClique) {
    22442205      // initialize new clique
    2245       abort();
     2206      lastRow=iRow;
     2207      cliquePointer=choleskyStart_[iRow];
     2208      cliqueSize=choleskyStart_[iRow+1]-cliquePointer+1;
    22462209    }
    22472210    // for each column L[*,kRow] that affects L[*,iRow]
    2248     double diagonalValue=diagonal_[iRow];
     2211    longDouble diagonalValue=diagonal_[iRow];
    22492212    int nextRow = link_[iRow];
    22502213    int kRow=0;
     
    22582221      CoinBigIndex end=choleskyStart_[kRow+1];
    22592222      assert(k<end);
    2260       double a_ik=sparseFactor_[k++];
    2261       double value1=d[kRow]*a_ik;
     2223      longDouble a_ik=sparseFactor_[k++];
     2224      longDouble value1=d[kRow]*a_ik;
    22622225      // update first
    22632226      first[kRow]=k;
    22642227      diagonalValue -= value1*a_ik;
     2228      CoinBigIndex offset = indexStart_[kRow]-choleskyStart_[kRow];
     2229      int jRow = choleskyRow_[k+offset];
    22652230      if (k<end) {
    2266         int offset = indexStart_[kRow]-choleskyStart_[kRow];
    2267         int jRow = choleskyRow_[k+offset];
    2268         link_[kRow]=link_[jRow];
    2269         link_[jRow]=kRow;
    22702231        if (clique_[kRow]<MINCLIQUE) {
     2232          link_[kRow]=link_[jRow];
     2233          link_[jRow]=kRow;
    22712234          for (;k<end;k++) {
    22722235            int jRow = choleskyRow_[k+offset];
    2273             work[jRow] -= sparseFactor_[k]*value1;
     2236            work[jRow] += sparseFactor_[k]*value1;
    22742237          }
    22752238        } else {
    22762239          // Clique
    2277           abort();
     2240          CoinBigIndex currentIndex = k+offset;
     2241          int linkSave=link_[jRow];
     2242          link_[jRow]=kRow;
     2243          work[kRow]=value1; // ? or a_jk
     2244          int last = kRow+clique_[kRow];
     2245          for (int kkRow=kRow+1;kkRow<last;kkRow++) {
     2246            CoinBigIndex j=first[kkRow];
     2247            longDouble a = sparseFactor_[j];
     2248            longDouble dValue = d[kkRow]*a;
     2249            diagonalValue -= a*dValue;
     2250            work[kkRow]=dValue;
     2251            first[kkRow]++;
     2252            link_[kkRow-1]=kkRow;
     2253          }
     2254          nextRow = link_[last-1];
     2255          link_[last-1]=linkSave;
     2256          int length = end-k+1;
     2257          for (int i=0;i<length;i++) {
     2258            int lRow = choleskyRow_[currentIndex++];
     2259            longDouble t0 = work[lRow];
     2260            for (int kkRow=kRow;kkRow<last;kkRow++) {
     2261              CoinBigIndex j = first[kkRow]+i;
     2262              t0 += work[kkRow]*sparseFactor_[j];
     2263            }
     2264            work[lRow]=t0;
     2265          }
    22782266        }
    22792267      }
     
    22822270    if (inClique) {
    22832271      // in clique
    2284       abort();
     2272      diagonal_[iRow]=diagonalValue;
     2273      CoinBigIndex start=choleskyStart_[iRow];
     2274      CoinBigIndex end=choleskyStart_[iRow+1];
     2275      CoinBigIndex currentIndex = indexStart_[iRow];
     2276      nextRow2=-1;
     2277      CoinBigIndex get=start+clique_[iRow]-1;
     2278      if (get<end) {
     2279        nextRow2 = choleskyRow_[currentIndex+get-start];
     2280        first[iRow]=get;
     2281      }
     2282      for (CoinBigIndex j=start;j<end;j++) {
     2283        int kRow = choleskyRow_[currentIndex++];
     2284        sparseFactor_[j] -= work[kRow]; // times?
     2285        work[kRow]=0.0;
     2286      }
    22852287    } else {
    22862288      // not in clique
     
    23142316      }
    23152317      diagonal_[iRow]=diagonalValue;
    2316       int offset = indexStart_[iRow]-choleskyStart_[iRow];
     2318      CoinBigIndex offset = indexStart_[iRow]-choleskyStart_[iRow];
    23172319      CoinBigIndex start = choleskyStart_[iRow];
    23182320      CoinBigIndex end = choleskyStart_[iRow+1];
     
    23242326        for (CoinBigIndex j=start;j<end;j++) {
    23252327          int jRow = choleskyRow_[j+offset];
    2326           double value = sparseFactor_[j]+work[jRow];
     2328          longDouble value = sparseFactor_[j]-work[jRow];
    23272329          work[jRow]=0.0;
    23282330          sparseFactor_[j]= diagonalValue*value;
     
    23312333    }
    23322334  }
    2333  
     2335  if (firstDense_<numberRows_) {
     2336    // do dense
     2337    // update dense part
     2338    updateDense(d,work,first);
     2339    ClpCholeskyDense dense;
     2340    // just borrow space
     2341    int nDense = numberRows_-firstDense_;
     2342    if (doKKT_) {
     2343      for (iRow=firstDense_;iRow<numberRows_;iRow++) {
     2344        int originalRow=permute_[iRow];
     2345        if (originalRow>=firstPositive) {
     2346          firstPositive=iRow-firstDense_;
     2347          break;
     2348        }
     2349      }
     2350    }
     2351    dense.reserveSpace(this,nDense);
     2352    int * dropped = new int[nDense];
     2353    memset(dropped,0,nDense*sizeof(int));
     2354    dense.setDoubleParameter(3,largest);
     2355    dense.setDoubleParameter(4,smallest);
     2356    dense.setDoubleParameter(10,dropValue);
     2357    dense.setIntegerParameter(20,0);
     2358    dense.setIntegerParameter(34,firstPositive);
     2359    dense.factorizePart2(dropped);
     2360    largest=dense.getDoubleParameter(3);
     2361    smallest=dense.getDoubleParameter(4);
     2362    integerParameters_[20]+=dense.getIntegerParameter(20);
     2363    for (iRow=firstDense_;iRow<numberRows_;iRow++) {
     2364      int originalRow=permute_[iRow];
     2365      rowsDropped[originalRow]=dropped[iRow-firstDense_];
     2366    }
     2367    delete [] dropped;
     2368  }
    23342369  delete [] d;
     2370  doubleParameters_[3]=largest;
     2371  doubleParameters_[4]=smallest;
    23352372  return;
     2373}
     2374// Updates dense part (broken out for profiling)
     2375void ClpCholeskyBase::updateDense(longDouble * d, longDouble * work, int * first)
     2376{
     2377  for (int iRow=0;iRow<firstDense_;iRow++) {
     2378    CoinBigIndex start=first[iRow];
     2379    CoinBigIndex end=choleskyStart_[iRow+1];
     2380    if (start<end) {
     2381      CoinBigIndex offset = indexStart_[iRow]-choleskyStart_[iRow];
     2382      longDouble dValue=d[iRow];
     2383      for (CoinBigIndex k=start;k<end;k++) {
     2384        int kRow = choleskyRow_[k+offset];
     2385        assert(kRow>=firstDense_);
     2386        longDouble a_ik=sparseFactor_[k];
     2387        longDouble value1=dValue*a_ik;
     2388        diagonal_[kRow] -= value1*a_ik;
     2389        CoinBigIndex base = choleskyStart_[kRow]-kRow-1;
     2390        for (CoinBigIndex j=k+1;j<end;j++) {
     2391          int jRow=choleskyRow_[j+offset];
     2392          longDouble a_jk = sparseFactor_[j];
     2393          sparseFactor_[base+jRow] -= a_jk*value1;
     2394        }
     2395      }
     2396    }
     2397  }
    23362398}
    23372399/* Uses factorization to solve. */
     
    23492411    double * change = new double[numberDense];
    23502412    for (i=0;i<numberDense;i++) {
    2351       const double * a = denseColumn_+i*numberRows_;
    2352       double value =0.0;
     2413      const longDouble * a = denseColumn_+i*numberRows_;
     2414      longDouble value =0.0;
    23532415      for (int iRow=0;iRow<numberRows_;iRow++)
    23542416        value += a[iRow]*region[iRow];
     
    23582420    dense_->solve(change);
    23592421    for (i=0;i<numberDense;i++) {
    2360       const double * a = denseColumn_+i*numberRows_;
    2361       double value = change[i];
     2422      const longDouble * a = denseColumn_+i*numberRows_;
     2423      longDouble value = change[i];
    23622424      for (int iRow=0;iRow<numberRows_;iRow++)
    23632425        region[iRow] -= value*a[iRow];
     
    23742436ClpCholeskyBase::solve(double * region, int type)
    23752437{
     2438#ifdef CLP_DEBUG
     2439  double * regionX=NULL;
     2440  if (sizeof(longWork)!=sizeof(double)&&type==3) {
     2441    regionX=ClpCopyOfArray(region,numberRows_);
     2442  }
     2443#endif
     2444  longWork * work = (longWork *) workDouble_;
     2445  int i;
     2446  CoinBigIndex j;
     2447  for (i=0;i<numberRows_;i++) {
     2448    int iRow = permute_[i];
     2449    work[i] = region[iRow];
     2450  }
     2451  switch (type) {
     2452  case 1:
     2453    for (i=0;i<numberRows_;i++) {
     2454      longDouble value=work[i];
     2455      CoinBigIndex offset = indexStart_[i]-choleskyStart_[i];
     2456      for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) {
     2457        int iRow = choleskyRow_[j+offset];
     2458        work[iRow] -= sparseFactor_[j]*value;
     2459      }
     2460    }
     2461    for (i=0;i<numberRows_;i++) {
     2462      int iRow = permute_[i];
     2463      region[iRow]=work[i]*diagonal_[i];
     2464    }
     2465    break;
     2466  case 2:
     2467    for (i=numberRows_-1;i>=0;i--) {
     2468      CoinBigIndex offset = indexStart_[i]-choleskyStart_[i];
     2469      longDouble value=work[i]*diagonal_[i];
     2470      for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) {
     2471        int iRow = choleskyRow_[j+offset];
     2472        value -= sparseFactor_[j]*work[iRow];
     2473      }
     2474      work[i]=value;
     2475      int iRow = permute_[i];
     2476      region[iRow]=value;
     2477    }
     2478    break;
     2479  case 3:
     2480    for (i=0;i<firstDense_;i++) {
     2481      CoinBigIndex offset = indexStart_[i]-choleskyStart_[i];
     2482      longDouble value=work[i];
     2483      for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) {
     2484        int iRow = choleskyRow_[j+offset];
     2485        work[iRow] -= sparseFactor_[j]*value;
     2486      }
     2487    }
     2488    if (firstDense_<numberRows_) {
     2489      // do dense
     2490      ClpCholeskyDense dense;
     2491      // just borrow space
     2492      int nDense = numberRows_-firstDense_;
     2493      dense.reserveSpace(this,nDense);
     2494#if CLP_LONG_CHOLESKY!=1
     2495      dense.solveLong(work+firstDense_);
     2496#else
     2497      dense.solveLongWork(work+firstDense_);
     2498#endif
     2499      for (i=numberRows_-1;i>=firstDense_;i--) {
     2500        longDouble value=work[i];
     2501        int iRow = permute_[i];
     2502        region[iRow]=value;
     2503      }
     2504    }
     2505    for (i=firstDense_-1;i>=0;i--) {
     2506      CoinBigIndex offset = indexStart_[i]-choleskyStart_[i];
     2507      longDouble value=work[i]*diagonal_[i];
     2508      for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) {
     2509        int iRow = choleskyRow_[j+offset];
     2510        value -= sparseFactor_[j]*work[iRow];
     2511      }
     2512      work[i]=value;
     2513      int iRow = permute_[i];
     2514      region[iRow]=value;
     2515    }
     2516    break;
     2517  }
     2518#ifdef CLP_DEBUG
     2519  if (regionX) {
     2520    longDouble * work = workDouble_;
     2521    int i;
     2522    CoinBigIndex j;
     2523    double largestO=0.0;
     2524    for (i=0;i<numberRows_;i++) {
     2525      largestO = CoinMax(largestO,fabs(regionX[i]));
     2526    }
     2527    for (i=0;i<numberRows_;i++) {
     2528      int iRow = permute_[i];
     2529      work[i] = regionX[iRow];
     2530    }
     2531    for (i=0;i<firstDense_;i++) {
     2532      CoinBigIndex offset = indexStart_[i]-choleskyStart_[i];
     2533      longDouble value=work[i];
     2534      for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) {
     2535        int iRow = choleskyRow_[j+offset];
     2536        work[iRow] -= sparseFactor_[j]*value;
     2537      }
     2538    }
     2539    if (firstDense_<numberRows_) {
     2540      // do dense
     2541      ClpCholeskyDense dense;
     2542      // just borrow space
     2543      int nDense = numberRows_-firstDense_;
     2544      dense.reserveSpace(this,nDense);
     2545      dense.solveLong(work+firstDense_);
     2546      for (i=numberRows_-1;i>=firstDense_;i--) {
     2547        longDouble value=work[i];
     2548        int iRow = permute_[i];
     2549        regionX[iRow]=value;
     2550      }
     2551    }
     2552    for (i=firstDense_-1;i>=0;i--) {
     2553      CoinBigIndex offset = indexStart_[i]-choleskyStart_[i];
     2554      longDouble value=work[i]*diagonal_[i];
     2555      for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) {
     2556        int iRow = choleskyRow_[j+offset];
     2557        value -= sparseFactor_[j]*work[iRow];
     2558      }
     2559      work[i]=value;
     2560      int iRow = permute_[i];
     2561      regionX[iRow]=value;
     2562    }
     2563    double largest=0.0;
     2564    double largestV=0.0;
     2565    for (i=0;i<numberRows_;i++) {
     2566      largest = CoinMax(largest,fabs(region[i]-regionX[i]));
     2567      largestV = CoinMax(largestV,fabs(region[i]));
     2568    }
     2569    printf("largest difference %g, largest %g, largest original %g\n",
     2570           largest,largestV,largestO);
     2571    delete [] regionX;
     2572  }
     2573#endif
     2574}
     2575/* solve - 1 just first half, 2 just second half - 3 both.
     2576   If 1 and 2 then diagonal has sqrt of inverse otherwise inverse
     2577*/
     2578void
     2579ClpCholeskyBase::solveLong(longDouble * region, int type)
     2580{
    23762581  int i;
    23772582  CoinBigIndex j;
     
    23832588  case 1:
    23842589    for (i=0;i<numberRows_;i++) {
    2385       double value=workDouble_[i];
     2590      longDouble value=workDouble_[i];
     2591      CoinBigIndex offset = indexStart_[i]-choleskyStart_[i];
    23862592      for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) {
    2387         int iRow = choleskyRow_[j];
     2593        int iRow = choleskyRow_[j+offset];
    23882594        workDouble_[iRow] -= sparseFactor_[j]*value;
    23892595      }
     
    23962602  case 2:
    23972603    for (i=numberRows_-1;i>=0;i--) {
    2398       double value=workDouble_[i]*diagonal_[i];
     2604      CoinBigIndex offset = indexStart_[i]-choleskyStart_[i];
     2605      longDouble value=workDouble_[i]*diagonal_[i];
    23992606      for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) {
    2400         int iRow = choleskyRow_[j];
     2607        int iRow = choleskyRow_[j+offset];
    24012608        value -= sparseFactor_[j]*workDouble_[iRow];
    24022609      }
     
    24072614    break;
    24082615  case 3:
    2409     for (i=0;i<numberRows_;i++) {
    2410       double value=workDouble_[i];
     2616    for (i=0;i<firstDense_;i++) {
     2617      CoinBigIndex offset = indexStart_[i]-choleskyStart_[i];
     2618      longDouble value=workDouble_[i];
    24112619      for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) {
    2412         int iRow = choleskyRow_[j];
     2620        int iRow = choleskyRow_[j+offset];
    24132621        workDouble_[iRow] -= sparseFactor_[j]*value;
    24142622      }
    24152623    }
    2416     for (i=numberRows_-1;i>=0;i--) {
    2417       double value=workDouble_[i]*diagonal_[i];
     2624    if (firstDense_<numberRows_) {
     2625      // do dense
     2626      ClpCholeskyDense dense;
     2627      // just borrow space
     2628      int nDense = numberRows_-firstDense_;
     2629      dense.reserveSpace(this,nDense);
     2630      dense.solveLong(workDouble_+firstDense_);
     2631      for (i=numberRows_-1;i>=firstDense_;i--) {
     2632        longDouble value=workDouble_[i];
     2633        int iRow = permute_[i];
     2634        region[iRow]=value;
     2635      }
     2636    }
     2637    for (i=firstDense_-1;i>=0;i--) {
     2638      CoinBigIndex offset = indexStart_[i]-choleskyStart_[i];
     2639      longDouble value=workDouble_[i]*diagonal_[i];
    24182640      for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) {
    2419         int iRow = choleskyRow_[j];
     2641        int iRow = choleskyRow_[j+offset];
    24202642        value -= sparseFactor_[j]*workDouble_[iRow];
    24212643      }
  • trunk/ClpCholeskyDense.cpp

    r414 r423  
    8888  numberBlocks = numberBlocks + ((numberBlocks*(numberBlocks+1))/2);
    8989  sizeFactor_=numberBlocks*BLOCKSQ;
    90 #define CHOL_COMPARE
     90  //#define CHOL_COMPARE
    9191#ifdef CHOL_COMPARE 
    9292  sizeFactor_ += 95000;
    9393#endif
    9494  if (!factor) {
    95     sparseFactor_ = new double [sizeFactor_];
     95    sparseFactor_ = new longDouble [sizeFactor_];
    9696    rowsDropped_ = new char [numberRows_];
    9797    memset(rowsDropped_,0,numberRows_);
    98     workDouble_ = new double[numberRows_];
    99     diagonal_ = new double[numberRows_];
     98    workDouble_ = new longDouble[numberRows_];
     99    diagonal_ = new longDouble[numberRows_];
    100100  } else {
    101101    borrowSpace_=true;
     
    163163  CoinZeroN(sparseFactor_,sizeFactor_);
    164164  //perturbation
    165   double perturbation=model_->diagonalPerturbation()*model_->diagonalNorm();
     165  longDouble perturbation=model_->diagonalPerturbation()*model_->diagonalNorm();
    166166  perturbation=perturbation*perturbation;
    167167  if (perturbation>1.0) {
     
    175175  double largest=1.0;
    176176  double smallest=COIN_DBL_MAX;
    177   double delta2 = model_->delta(); // add delta*delta to diagonal
     177  longDouble delta2 = model_->delta(); // add delta*delta to diagonal
    178178  delta2 *= delta2;
    179179  if (!doKKT_) {
    180     double * work = sparseFactor_;
     180    longDouble * work = sparseFactor_;
    181181    work--; // skip diagonal
    182182    int addOffset=numberRows_-1;
     
    188188        CoinBigIndex startRow=rowStart[iRow];
    189189        CoinBigIndex endRow=rowStart[iRow]+rowLength[iRow];
    190         double diagonalValue = diagonalSlack[iRow]+delta2;
     190        longDouble diagonalValue = diagonalSlack[iRow]+delta2;
    191191        for (CoinBigIndex k=startRow;k<endRow;k++) {
    192192          int iColumn=column[k];
    193193          CoinBigIndex start=columnStart[iColumn];
    194194          CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
    195           double multiplier = diagonal[iColumn]*elementByRow[k];
     195          longDouble multiplier = diagonal[iColumn]*elementByRow[k];
    196196          for (CoinBigIndex j=start;j<end;j++) {
    197197            int jRow=row[j];
    198198            if (!rowsDropped_[jRow]) {
    199199              if (jRow>iRow) {
    200                 double value=element[j]*multiplier;
     200                longDouble value=element[j]*multiplier;
    201201                work[jRow] += value;
    202202              } else if (jRow==iRow) {
    203                 double value=element[j]*multiplier;
     203                longDouble value=element[j]*multiplier;
    204204                diagonalValue += value;
    205205              }
     
    227227      rowsDropped[iRow]=dropped;
    228228      if (!dropped) {
    229         double diagonal = diagonal_[iRow];
     229        longDouble diagonal = diagonal_[iRow];
    230230        if (diagonal>largest2) {
    231231          diagonal_[iRow]=diagonal+perturbation;
     
    279279    int numberColumns = model_->numberColumns();
    280280    int numberTotal = numberColumns + numberRowsModel;
    281     double * work = sparseFactor_;
     281    longDouble * work = sparseFactor_;
    282282    work--; // skip diagonal
    283283    int addOffset=numberRows_-1;
     
    285285    if (!quadratic) {
    286286      for (iColumn=0;iColumn<numberColumns;iColumn++) {
    287         double value = diagonal[iColumn];
     287        longDouble value = diagonal[iColumn];
    288288        if (fabs(value)>1.0e-100) {
    289289          value = 1.0/value;
     
    311311      const double * quadraticElement = quadratic->getElements();
    312312      for (iColumn=0;iColumn<numberColumns;iColumn++) {
    313         double value = diagonal[iColumn];
     313        longDouble value = diagonal[iColumn];
    314314        CoinBigIndex j;
    315315        if (fabs(value)>1.0e-100) {
     
    342342    // slacks
    343343    for (iColumn=numberColumns;iColumn<numberTotal;iColumn++) {
    344       double value = diagonal[iColumn];
     344      longDouble value = diagonal[iColumn];
    345345      if (fabs(value)>1.0e-100) {
    346346        value = 1.0/value;
     
    417417{
    418418  int iColumn;
    419   double * xx = sparseFactor_;
    420   double * yy = diagonal_;
     419  longDouble * xx = sparseFactor_;
     420  longDouble * yy = diagonal_;
    421421  diagonal_ = sparseFactor_ + 40000;
    422422  sparseFactor_=diagonal_ + numberRows_;
     
    429429  double dropValue = doubleParameters_[10];
    430430  int firstPositive=integerParameters_[34];
    431   double * work = sparseFactor_;
     431  longDouble * work = sparseFactor_;
    432432  // Allow for triangular
    433433  int addOffset=numberRows_-1;
     
    436436    int iRow;
    437437    int addOffsetNow = numberRows_-1;;
    438     double * workNow=sparseFactor_-1+iColumn;
     438    longDouble * workNow=sparseFactor_-1+iColumn;
    439439    double diagonalValue = diagonal_[iColumn];
    440440    for (iRow=0;iRow<iColumn;iRow++) {
     
    510510#ifdef POS_DEBUG
    511511static int counter=0;
    512 int ClpCholeskyDense::bNumber(const double * array, int &iRow, int &iCol)
     512int ClpCholeskyDense::bNumber(const longDouble * array, int &iRow, int &iCol)
    513513{
    514514  int numberBlocks = (numberRows_+BLOCK-1)>>BLOCKSHIFT;
    515   double * a = sparseFactor_+BLOCKSQ*numberBlocks;
     515  longDouble * a = sparseFactor_+BLOCKSQ*numberBlocks;
    516516  int k = array-a;
    517517  assert ((k%BLOCKSQ)==0);
     
    537537{
    538538  int iColumn;
    539   //double * xx = sparseFactor_;
    540   //double * yy = diagonal_;
     539  //longDouble * xx = sparseFactor_;
     540  //longDouble * yy = diagonal_;
    541541  //diagonal_ = sparseFactor_ + 40000;
    542542  //sparseFactor_=diagonal_ + numberRows_;
     
    546546  int numberBlocks = (numberRows_+BLOCK-1)>>BLOCKSHIFT;
    547547  // later align on boundary
    548   double * a = sparseFactor_+BLOCKSQ*numberBlocks;
     548  longDouble * a = sparseFactor_+BLOCKSQ*numberBlocks;
    549549  int n=numberRows_;
    550550  int nRound = numberRows_&(~(BLOCK-1));
     
    558558  int rowLast;
    559559  if (sizeLastBlock!=BLOCK) {
    560     double * aa = &a[(block-1)*BLOCKSQ];
     560    longDouble * aa = &a[(block-1)*BLOCKSQ];
    561561    rowLast=nRound-1;
    562562    ifOdd=1;
     
    584584  int nBlock=0;
    585585  for (;n>0;n-=BLOCK) {
    586     double * aa = &a[(block-1)*BLOCKSQ];
    587     double * aaLast=NULL;
     586    longDouble * aa = &a[(block-1)*BLOCKSQ];
     587    longDouble * aaLast=NULL;
    588588    int put = BLOCKSQ;
    589589    int putLast=0;
     
    603603        putLast -= BLOCK-sizeLastBlock;
    604604      }
    605       double * aPut = aa;
     605      longDouble * aPut = aa;
    606606      int j=rowLast;
    607607      for (int jBlock=0;jBlock<=nBlock;jBlock++) {
     
    631631// Non leaf recursive factor
    632632void
    633 ClpCholeskyDense::factor(double * a, int n, int numberBlocks,
    634             double * diagonal, double * work, int * rowsDropped)
     633ClpCholeskyDense::factor(longDouble * a, int n, int numberBlocks,
     634            longDouble * diagonal, longDouble * work, int * rowsDropped)
    635635{
    636636  if (n<=BLOCK) {
     
    639639    int nb=number_blocks((n+1)>>1);
    640640    int nThis=number_rows(nb);
    641     double * aother;
     641    longDouble * aother;
    642642    int nLeft=n-nThis;
    643643    int nintri=(nb*(nb+1))>>1;
     
    653653// Non leaf recursive triangle rectangle update
    654654void
    655 ClpCholeskyDense::triRec(double * aTri, int nThis, double * aUnder,
    656                          double * diagonal, double * work,
     655ClpCholeskyDense::triRec(longDouble * aTri, int nThis, longDouble * aUnder,
     656                         longDouble * diagonal, longDouble * work,
    657657                         int nLeft, int iBlock, int jBlock,
    658658                         int numberBlocks)
     
    669669    int nb=number_blocks((nThis+1)>>1);
    670670    int nThis2=number_rows(nb);
    671     double * aother;
     671    longDouble * aother;
    672672    int kBlock=jBlock+nb;
    673673    int i;
     
    688688// Non leaf recursive rectangle triangle update
    689689void
    690 ClpCholeskyDense::recTri(double * aUnder, int nTri, int nDo,
    691                          int iBlock, int jBlock,double * aTri,
    692                          double * diagonal, double * work,
     690ClpCholeskyDense::recTri(longDouble * aUnder, int nTri, int nDo,
     691                         int iBlock, int jBlock,longDouble * aTri,
     692                         longDouble * diagonal, longDouble * work,
    693693                         int numberBlocks)
    694694{
     
    698698    int nb=number_blocks((nDo+1)>>1);
    699699    int nDo2=number_rows(nb);
    700     double * aother;
     700    longDouble * aother;
    701701    int i;
    702702    recTri(aUnder,nTri,nDo2,iBlock,jBlock,aTri,diagonal,work,numberBlocks);
     
    709709    int nb=number_blocks((nTri+1)>>1);
    710710    int nTri2=number_rows(nb);
    711     double * aother;
     711    longDouble * aother;
    712712    int i;
    713713    recTri(aUnder,nTri2,nDo,iBlock,jBlock,aTri,diagonal,work,numberBlocks);
     
    727727*/
    728728void
    729 ClpCholeskyDense::recRec(double * above, int nUnder, int nUnderK,
    730                          int nDo, double * aUnder, double *aOther, double * diagonal, double * work,
     729ClpCholeskyDense::recRec(longDouble * above, int nUnder, int nUnderK,
     730                         int nDo, longDouble * aUnder, longDouble *aOther, longDouble * diagonal, longDouble * work,
    731731            int kBlock,int iBlock, int jBlock,
    732732            int numberBlocks)
     
    767767// Leaf recursive factor
    768768void
    769 ClpCholeskyDense::factorLeaf(double * a, int n,
    770                 double * diagonal, double * work, int * rowsDropped)
    771 {
    772   double largest=doubleParameters_[3];
    773   double smallest=doubleParameters_[4];
     769ClpCholeskyDense::factorLeaf(longDouble * a, int n,
     770                longDouble * diagonal, longDouble * work, int * rowsDropped)
     771{
     772  longDouble largest=doubleParameters_[3];
     773  longDouble smallest=doubleParameters_[4];
    774774  double dropValue = doubleParameters_[10];
    775775  int firstPositive=integerParameters_[34];
     
    777777  int numberDropped=0;
    778778  int i, j, k;
    779   double t00,temp1, * aa;
     779  longDouble t00,temp1;
     780  longDouble * aa;
    780781  aa = a-BLOCK;
    781782  for (j = 0; j < n; j ++) {
     
    783784    t00 = aa[j];
    784785    for (k = 0; k < j; ++k) {
    785       double multiplier = work[k];
     786      longDouble multiplier = work[k];
    786787      t00 -= a[j + k * BLOCK] * a[j + k * BLOCK] * multiplier;
    787788    }
    788789    bool dropColumn=false;
    789     double useT00=t00;
     790    longDouble useT00=t00;
    790791    if (j+rowOffset<firstPositive) {
    791792      // must be negative
     
    824825        t00=aa[i];
    825826        for (k = 0; k < j; ++k) {
    826           double multiplier = work[k];
     827          longDouble multiplier = work[k];
    827828          t00 -= a[i + k * BLOCK] * a[j + k * BLOCK] * multiplier;
    828829        }
     
    847848// Leaf recursive triangle rectangle update
    848849void
    849 ClpCholeskyDense::triRecLeaf(double * aTri, double * aUnder, double * diagonal, double * work,
     850ClpCholeskyDense::triRecLeaf(longDouble * aTri, longDouble * aUnder, longDouble * diagonal, longDouble * work,
    850851                int nUnder)
    851852{
     
    861862#endif
    862863  int j;
    863   double * aa;
     864  longDouble * aa;
    864865#if BLOCK==16
    865866  if (nUnder==BLOCK) {
     
    867868    for (j = 0; j < BLOCK; j +=2) {
    868869      aa+=2*BLOCK;
    869       double temp0 = diagonal[j];
    870       double temp1 = diagonal[j+1];
     870      longDouble temp0 = diagonal[j];
     871      longDouble temp1 = diagonal[j+1];
    871872      for (int i=0;i<BLOCK;i+=2) {
    872         double t00=aUnder[i+j*BLOCK];
    873         double t10=aUnder[i+ BLOCK + j*BLOCK];
    874         double t01=aUnder[i+1+j*BLOCK];
    875         double t11=aUnder[i+1+ BLOCK + j*BLOCK];
     873        longDouble t00=aUnder[i+j*BLOCK];
     874        longDouble t10=aUnder[i+ BLOCK + j*BLOCK];
     875        longDouble t01=aUnder[i+1+j*BLOCK];
     876        longDouble t11=aUnder[i+1+ BLOCK + j*BLOCK];
    876877        for (int k = 0; k < j; ++k) {
    877           double multiplier=work[k];
    878           double au0 = aUnder[i + k * BLOCK]*multiplier;
    879           double au1 = aUnder[i + 1 + k * BLOCK]*multiplier;
    880           double at0 = aTri[j + k * BLOCK];
    881           double at1 = aTri[j + 1 + k * BLOCK];
     878          longDouble multiplier=work[k];
     879          longDouble au0 = aUnder[i + k * BLOCK]*multiplier;
     880          longDouble au1 = aUnder[i + 1 + k * BLOCK]*multiplier;
     881          longDouble at0 = aTri[j + k * BLOCK];
     882          longDouble at1 = aTri[j + 1 + k * BLOCK];
    882883          t00 -= au0 * at0;
    883884          t10 -= au0 * at1;
     
    886887        }
    887888        t00 *= temp0;
    888         double at1 = aTri[j + 1 + j*BLOCK]*work[j];
     889        longDouble at1 = aTri[j + 1 + j*BLOCK]*work[j];
    889890        t10 -= t00 * at1;
    890891        t01 *= temp0;
     
    901902    for (j = 0; j < BLOCK; j ++) {
    902903      aa+=BLOCK;
    903       double temp1 = diagonal[j];
     904      longDouble temp1 = diagonal[j];
    904905      for (int i=0;i<nUnder;i++) {
    905         double t00=aUnder[i+j*BLOCK];
     906        longDouble t00=aUnder[i+j*BLOCK];
    906907        for (int k = 0; k < j; ++k) {
    907           double multiplier=work[k];
     908          longDouble multiplier=work[k];
    908909          t00 -= aUnder[i + k * BLOCK] * aTri[j + k * BLOCK] * multiplier;
    909910        }
     
    916917}
    917918// Leaf recursive rectangle triangle update
    918 void ClpCholeskyDense::recTriLeaf(double * aUnder, double * aTri,
    919                                   double * diagonal, double * work, int nUnder)
     919void ClpCholeskyDense::recTriLeaf(longDouble * aUnder, longDouble * aTri,
     920                                  longDouble * diagonal, longDouble * work, int nUnder)
    920921{
    921922#ifdef POS_DEBUG
     
    930931#endif
    931932  int i, j, k;
    932   double t00, * aa;
     933  longDouble t00;
     934  longDouble * aa;
    933935#if BLOCK==16
    934936  if (nUnder==BLOCK) {
    935     double * aUnder2 = aUnder-2;
     937    longDouble * aUnder2 = aUnder-2;
    936938    aa = aTri-2*BLOCK;
    937939    for (j = 0; j < BLOCK; j +=2) {
    938       double t00,t01,t10,t11;
     940      longDouble t00,t01,t10,t11;
    939941      aa+=2*BLOCK;
    940942      aUnder2+=2;
     
    943945      t10=aa[j+1+BLOCK];
    944946      for (k = 0; k < BLOCK; ++k) {
    945         double multiplier = work[k];
    946         double a0 = aUnder2[k * BLOCK];
    947         double a1 = aUnder2[1 + k * BLOCK];
    948         double x0 = a0 * multiplier;
    949         double x1 = a1 * multiplier;
     947        longDouble multiplier = work[k];
     948        longDouble a0 = aUnder2[k * BLOCK];
     949        longDouble a1 = aUnder2[1 + k * BLOCK];
     950        longDouble x0 = a0 * multiplier;
     951        longDouble x1 = a1 * multiplier;
    950952        t00 -= a0 * x0;
    951953        t01 -= a1 * x0;
     
    961963        t11=aa[i+1+BLOCK];
    962964        for (k = 0; k < BLOCK; ++k) {
    963           double multiplier = work[k];
    964           double a0 = aUnder2[k * BLOCK]*multiplier;
    965           double a1 = aUnder2[1 + k * BLOCK]*multiplier;
     965          longDouble multiplier = work[k];
     966          longDouble a0 = aUnder2[k * BLOCK]*multiplier;
     967          longDouble a1 = aUnder2[1 + k * BLOCK]*multiplier;
    966968          t00 -= aUnder[i + k * BLOCK] * a0;
    967969          t01 -= aUnder[i + k * BLOCK] * a1;
     
    983985        t00=aa[i];
    984986        for (k = 0; k < BLOCK; ++k) {
    985           double multiplier = work[k];
     987          longDouble multiplier = work[k];
    986988          t00 -= aUnder[i + k * BLOCK] * aUnder[j + k * BLOCK]*multiplier;
    987989        }
     
    9981000*/
    9991001void
    1000 ClpCholeskyDense::recRecLeaf(double * above,
    1001                              double * aUnder, double *aOther, double * diagonal, double * work,
     1002ClpCholeskyDense::recRecLeaf(longDouble * above,
     1003                             longDouble * aUnder, longDouble *aOther, longDouble * diagonal, longDouble * work,
    10021004                             int nUnder)
    10031005{
     
    10151017#endif
    10161018  int i, j, k;
    1017   double * aa;
     1019  longDouble * aa;
    10181020#if BLOCK==16
    10191021  aa = aOther-4*BLOCK;
     
    10221024      aa+=4*BLOCK;
    10231025      for (i=0;i< BLOCK;i+=4) {
    1024         double t00=aa[i+0+0*BLOCK];
    1025         double t10=aa[i+0+1*BLOCK];
    1026         double t20=aa[i+0+2*BLOCK];
    1027         double t30=aa[i+0+3*BLOCK];
    1028         double t01=aa[i+1+0*BLOCK];
    1029         double t11=aa[i+1+1*BLOCK];
    1030         double t21=aa[i+1+2*BLOCK];
    1031         double t31=aa[i+1+3*BLOCK];
    1032         double t02=aa[i+2+0*BLOCK];
    1033         double t12=aa[i+2+1*BLOCK];
    1034         double t22=aa[i+2+2*BLOCK];
    1035         double t32=aa[i+2+3*BLOCK];
    1036         double t03=aa[i+3+0*BLOCK];
    1037         double t13=aa[i+3+1*BLOCK];
    1038         double t23=aa[i+3+2*BLOCK];
    1039         double t33=aa[i+3+3*BLOCK];
     1026        longDouble t00=aa[i+0+0*BLOCK];
     1027        longDouble t10=aa[i+0+1*BLOCK];
     1028        longDouble t20=aa[i+0+2*BLOCK];
     1029        longDouble t30=aa[i+0+3*BLOCK];
     1030        longDouble t01=aa[i+1+0*BLOCK];
     1031        longDouble t11=aa[i+1+1*BLOCK];
     1032        longDouble t21=aa[i+1+2*BLOCK];
     1033        longDouble t31=aa[i+1+3*BLOCK];
     1034        longDouble t02=aa[i+2+0*BLOCK];
     1035        longDouble t12=aa[i+2+1*BLOCK];
     1036        longDouble t22=aa[i+2+2*BLOCK];
     1037        longDouble t32=aa[i+2+3*BLOCK];
     1038        longDouble t03=aa[i+3+0*BLOCK];
     1039        longDouble t13=aa[i+3+1*BLOCK];
     1040        longDouble t23=aa[i+3+2*BLOCK];
     1041        longDouble t33=aa[i+3+3*BLOCK];
    10401042        for (k=0;k<BLOCK;k++) {
    1041           double multiplier = work[k];
    1042           double a00=aUnder[i+0+k*BLOCK]*multiplier;
    1043           double a01=aUnder[i+1+k*BLOCK]*multiplier;
    1044           double a02=aUnder[i+2+k*BLOCK]*multiplier;
    1045           double a03=aUnder[i+3+k*BLOCK]*multiplier;
     1043          longDouble multiplier = work[k];
     1044          longDouble a00=aUnder[i+0+k*BLOCK]*multiplier;
     1045          longDouble a01=aUnder[i+1+k*BLOCK]*multiplier;
     1046          longDouble a02=aUnder[i+2+k*BLOCK]*multiplier;
     1047          longDouble a03=aUnder[i+3+k*BLOCK]*multiplier;
    10461048          t00 -= a00 * above[j + 0 + k * BLOCK];
    10471049          t10 -= a00 * above[j + 1 + k * BLOCK];
     
    10851087      aa+=4*BLOCK;
    10861088      for (i=0;i< n;i+=2) {
    1087         double t00=aa[i+0*BLOCK];
    1088         double t10=aa[i+1*BLOCK];
    1089         double t20=aa[i+2*BLOCK];
    1090         double t30=aa[i+3*BLOCK];
    1091         double t01=aa[i+1+0*BLOCK];
    1092         double t11=aa[i+1+1*BLOCK];
    1093         double t21=aa[i+1+2*BLOCK];
    1094         double t31=aa[i+1+3*BLOCK];
     1089        longDouble t00=aa[i+0*BLOCK];
     1090        longDouble t10=aa[i+1*BLOCK];
     1091        longDouble t20=aa[i+2*BLOCK];
     1092        longDouble t30=aa[i+3*BLOCK];
     1093        longDouble t01=aa[i+1+0*BLOCK];
     1094        longDouble t11=aa[i+1+1*BLOCK];
     1095        longDouble t21=aa[i+1+2*BLOCK];
     1096        longDouble t31=aa[i+1+3*BLOCK];
    10951097        for (k=0;k<BLOCK;k++) {
    1096           double multiplier = work[k];
    1097           double a00=aUnder[i+k*BLOCK]*multiplier;
    1098           double a01=aUnder[i+1+k*BLOCK]*multiplier;
     1098          longDouble multiplier = work[k];
     1099          longDouble a00=aUnder[i+k*BLOCK]*multiplier;
     1100          longDouble a01=aUnder[i+1+k*BLOCK]*multiplier;
    10991101          t00 -= a00 * above[j + 0 + k * BLOCK];
    11001102          t10 -= a00 * above[j + 1 + k * BLOCK];
     
    11161118      }
    11171119      if (odd) {
    1118         double t0=aa[n+0*BLOCK];
    1119         double t1=aa[n+1*BLOCK];
    1120         double t2=aa[n+2*BLOCK];
    1121         double t3=aa[n+3*BLOCK];
    1122         double a0;
     1120        longDouble t0=aa[n+0*BLOCK];
     1121        longDouble t1=aa[n+1*BLOCK];
     1122        longDouble t2=aa[n+2*BLOCK];
     1123        longDouble t3=aa[n+3*BLOCK];
     1124        longDouble a0;
    11231125        for (k=0;k<BLOCK;k++) {
    11241126          a0=aUnder[n+k*BLOCK]*work[k];
     
    11401142    aa+=BLOCK;
    11411143    for (i=0;i< nUnder;i++) {
    1142       double t00=aa[i+0*BLOCK];
     1144      longDouble t00=aa[i+0*BLOCK];
    11431145      for (k=0;k<BLOCK;k++) {
    1144         double a00=aUnder[i+k*BLOCK]*work[k];
     1146        longDouble a00=aUnder[i+k*BLOCK]*work[k];
    11451147        t00 -= a00 * above[j + k * BLOCK];
    11461148      }
     
    11571159  double * region2=NULL;
    11581160  if (numberRows_<200) {
    1159     double * xx = sparseFactor_;
    1160     double * yy = diagonal_;
     1161    longDouble * xx = sparseFactor_;
     1162    longDouble * yy = diagonal_;
    11611163    diagonal_ = sparseFactor_ + 40000;
    11621164    sparseFactor_=diagonal_ + numberRows_;
     
    11641166    int iRow,iColumn;
    11651167    int addOffset=numberRows_-1;
    1166     double * work = sparseFactor_-1;
     1168    longDouble * work = sparseFactor_-1;
    11671169    for (iColumn=0;iColumn<numberRows_;iColumn++) {
    11681170      double value = region2[iColumn];
     
    11841186  }
    11851187#endif
    1186   //double * xx = sparseFactor_;
    1187   //double * yy = diagonal_;
     1188  //longDouble * xx = sparseFactor_;
     1189  //longDouble * yy = diagonal_;
    11881190  //diagonal_ = sparseFactor_ + 40000;
    11891191  //sparseFactor_=diagonal_ + numberRows_;
    11901192  int numberBlocks = (numberRows_+BLOCK-1)>>BLOCKSHIFT;
    11911193  // later align on boundary
    1192   double * a = sparseFactor_+BLOCKSQ*numberBlocks;
     1194  longDouble * a = sparseFactor_+BLOCKSQ*numberBlocks;
    11931195  int iBlock;
    1194   double * aa = a;
     1196  longDouble * aa = a;
    11951197  int iColumn;
    11961198  for (iBlock=0;iBlock<numberBlocks;iBlock++) {
     
    12571259// Forward part of solve 1
    12581260void
    1259 ClpCholeskyDense::solveF1(double * a,int n,double * region)
     1261ClpCholeskyDense::solveF1(longDouble * a,int n,double * region)
    12601262{
    12611263  int j, k;
    1262   double t00;
     1264  longDouble t00;
    12631265  for (j = 0; j < n; j ++) {
    12641266    t00 = region[j];
     
    12721274// Forward part of solve 2
    12731275void
    1274 ClpCholeskyDense::solveF2(double * a,int n,double * region, double * region2)
     1276ClpCholeskyDense::solveF2(longDouble * a,int n,double * region, double * region2)
    12751277{
    12761278  int j, k;
     
    12781280  if (n==BLOCK) {
    12791281    for (k = 0; k < BLOCK; k+=4) {
    1280       double t0 = region2[0];
    1281       double t1 = region2[1];
    1282       double t2 = region2[2];
    1283       double t3 = region2[3];
     1282      longDouble t0 = region2[0];
     1283      longDouble t1 = region2[1];
     1284      longDouble t2 = region2[2];
     1285      longDouble t3 = region2[3];
    12841286      t0 -= region[0] * a[0 + 0 * BLOCK];
    12851287      t1 -= region[0] * a[1 + 0 * BLOCK];
     
    13721374#endif
    13731375    for (k = 0; k < n; ++k) {
    1374       double t00 = region2[k];
     1376      longDouble t00 = region2[k];
    13751377      for (j = 0; j < BLOCK; j ++) {
    13761378        t00 -= region[j] * a[k + j * BLOCK];
     
    13841386// Backward part of solve 1
    13851387void
    1386 ClpCholeskyDense::solveB1(double * a,int n,double * region)
     1388ClpCholeskyDense::solveB1(longDouble * a,int n,double * region)
    13871389{
    13881390  int j, k;
    1389   double t00;
     1391  longDouble t00;
    13901392  for (j = n-1; j >=0; j --) {
    13911393    t00 = region[j];
     
    13991401// Backward part of solve 2
    14001402void
    1401 ClpCholeskyDense::solveB2(double * a,int n,double * region, double * region2)
     1403ClpCholeskyDense::solveB2(longDouble * a,int n,double * region, double * region2)
    14021404{
    14031405  int j, k;
     
    14051407  if (n==BLOCK) {
    14061408    for (j = 0; j < BLOCK; j +=4) {
    1407       double t0 = region[0];
    1408       double t1 = region[1];
    1409       double t2 = region[2];
    1410       double t3 = region[3];
     1409      longDouble t0 = region[0];
     1410      longDouble t1 = region[1];
     1411      longDouble t2 = region[2];
     1412      longDouble t3 = region[3];
    14111413      t0 -= region2[0] * a[0 + 0*BLOCK];
    14121414      t1 -= region2[0] * a[0 + 1*BLOCK];
     
    15001502#endif
    15011503    for (j = 0; j < BLOCK; j ++) {
    1502       double t00 = region[j];
     1504      longDouble t00 = region[j];
    15031505      for (k = 0; k < n; ++k) {
    15041506        t00 -= region2[k] * a[k + j * BLOCK];
     
    15101512#endif
    15111513}
     1514/* Uses factorization to solve. */
     1515void
     1516ClpCholeskyDense::solveLong (longDouble * region)
     1517{
     1518  int numberBlocks = (numberRows_+BLOCK-1)>>BLOCKSHIFT;
     1519  // later align on boundary
     1520  longDouble * a = sparseFactor_+BLOCKSQ*numberBlocks;
     1521  int iBlock;
     1522  longDouble * aa = a;
     1523  int iColumn;
     1524  for (iBlock=0;iBlock<numberBlocks;iBlock++) {
     1525    int nChunk;
     1526    int jBlock;
     1527    int iDo = iBlock*BLOCK;
     1528    int base=iDo;
     1529    if (iDo+BLOCK>numberRows_) {
     1530      nChunk=numberRows_-iDo;
     1531    } else {
     1532      nChunk=BLOCK;
     1533    }
     1534    solveF1Long(aa,nChunk,region+iDo);
     1535    for (jBlock=iBlock+1;jBlock<numberBlocks;jBlock++) {
     1536      base+=BLOCK;
     1537      aa+=BLOCKSQ;
     1538      if (base+BLOCK>numberRows_) {
     1539        nChunk=numberRows_-base;
     1540      } else {
     1541        nChunk=BLOCK;
     1542      }
     1543      solveF2Long(aa,nChunk,region+iDo,region+base);
     1544    }
     1545    aa+=BLOCKSQ;
     1546  }
     1547  // do diagonal outside
     1548  for (iColumn=0;iColumn<numberRows_;iColumn++)
     1549    region[iColumn] *= diagonal_[iColumn];
     1550  int offset=((numberBlocks*(numberBlocks+1))>>1);
     1551  aa = a+number_entries(offset-1);
     1552  int lBase=(numberBlocks-1)*BLOCK;
     1553  for (iBlock=numberBlocks-1;iBlock>=0;iBlock--) {
     1554    int nChunk;
     1555    int jBlock;
     1556    int triBase=iBlock*BLOCK;
     1557    int iBase=lBase;
     1558    for (jBlock=iBlock+1;jBlock<numberBlocks;jBlock++) {
     1559      if (iBase+BLOCK>numberRows_) {
     1560        nChunk=numberRows_-iBase;
     1561      } else {
     1562        nChunk=BLOCK;
     1563      }
     1564      solveB2Long(aa,nChunk,region+triBase,region+iBase);
     1565      iBase-=BLOCK;
     1566      aa-=BLOCKSQ;
     1567    }
     1568    if (triBase+BLOCK>numberRows_) {
     1569      nChunk=numberRows_-triBase;
     1570    } else {
     1571      nChunk=BLOCK;
     1572    }
     1573    solveB1Long(aa,nChunk,region+triBase);
     1574    aa-=BLOCKSQ;
     1575  }
     1576}
     1577// Forward part of solve 1
     1578void
     1579ClpCholeskyDense::solveF1Long(longDouble * a,int n,longDouble * region)
     1580{
     1581  int j, k;
     1582  longDouble t00;
     1583  for (j = 0; j < n; j ++) {
     1584    t00 = region[j];
     1585    for (k = 0; k < j; ++k) {
     1586      t00 -= region[k] * a[j + k * BLOCK];
     1587    }
     1588    //t00*=a[j + j * BLOCK];
     1589    region[j] = t00;
     1590  }
     1591}
     1592// Forward part of solve 2
     1593void
     1594ClpCholeskyDense::solveF2Long(longDouble * a,int n,longDouble * region, longDouble * region2)
     1595{
     1596  int j, k;
     1597#if BLOCK==16
     1598  if (n==BLOCK) {
     1599    for (k = 0; k < BLOCK; k+=4) {
     1600      longDouble t0 = region2[0];
     1601      longDouble t1 = region2[1];
     1602      longDouble t2 = region2[2];
     1603      longDouble t3 = region2[3];
     1604      t0 -= region[0] * a[0 + 0 * BLOCK];
     1605      t1 -= region[0] * a[1 + 0 * BLOCK];
     1606      t2 -= region[0] * a[2 + 0 * BLOCK];
     1607      t3 -= region[0] * a[3 + 0 * BLOCK];
     1608
     1609      t0 -= region[1] * a[0 + 1 * BLOCK];
     1610      t1 -= region[1] * a[1 + 1 * BLOCK];
     1611      t2 -= region[1] * a[2 + 1 * BLOCK];
     1612      t3 -= region[1] * a[3 + 1 * BLOCK];
     1613
     1614      t0 -= region[2] * a[0 + 2 * BLOCK];
     1615      t1 -= region[2] * a[1 + 2 * BLOCK];
     1616      t2 -= region[2] * a[2 + 2 * BLOCK];
     1617      t3 -= region[2] * a[3 + 2 * BLOCK];
     1618
     1619      t0 -= region[3] * a[0 + 3 * BLOCK];
     1620      t1 -= region[3] * a[1 + 3 * BLOCK];
     1621      t2 -= region[3] * a[2 + 3 * BLOCK];
     1622      t3 -= region[3] * a[3 + 3 * BLOCK];
     1623
     1624      t0 -= region[4] * a[0 + 4 * BLOCK];
     1625      t1 -= region[4] * a[1 + 4 * BLOCK];
     1626      t2 -= region[4] * a[2 + 4 * BLOCK];
     1627      t3 -= region[4] * a[3 + 4 * BLOCK];
     1628
     1629      t0 -= region[5] * a[0 + 5 * BLOCK];
     1630      t1 -= region[5] * a[1 + 5 * BLOCK];
     1631      t2 -= region[5] * a[2 + 5 * BLOCK];
     1632      t3 -= region[5] * a[3 + 5 * BLOCK];
     1633
     1634      t0 -= region[6] * a[0 + 6 * BLOCK];
     1635      t1 -= region[6] * a[1 + 6 * BLOCK];
     1636      t2 -= region[6] * a[2 + 6 * BLOCK];
     1637      t3 -= region[6] * a[3 + 6 * BLOCK];
     1638
     1639      t0 -= region[7] * a[0 + 7 * BLOCK];
     1640      t1 -= region[7] * a[1 + 7 * BLOCK];
     1641      t2 -= region[7] * a[2 + 7 * BLOCK];
     1642      t3 -= region[7] * a[3 + 7 * BLOCK];
     1643#if BLOCK>8
     1644      t0 -= region[8] * a[0 + 8 * BLOCK];
     1645      t1 -= region[8] * a[1 + 8 * BLOCK];
     1646      t2 -= region[8] * a[2 + 8 * BLOCK];
     1647      t3 -= region[8] * a[3 + 8 * BLOCK];
     1648
     1649      t0 -= region[9] * a[0 + 9 * BLOCK];
     1650      t1 -= region[9] * a[1 + 9 * BLOCK];
     1651      t2 -= region[9] * a[2 + 9 * BLOCK];
     1652      t3 -= region[9] * a[3 + 9 * BLOCK];
     1653
     1654      t0 -= region[10] * a[0 + 10 * BLOCK];
     1655      t1 -= region[10] * a[1 + 10 * BLOCK];
     1656      t2 -= region[10] * a[2 + 10 * BLOCK];
     1657      t3 -= region[10] * a[3 + 10 * BLOCK];
     1658
     1659      t0 -= region[11] * a[0 + 11 * BLOCK];
     1660      t1 -= region[11] * a[1 + 11 * BLOCK];
     1661      t2 -= region[11] * a[2 + 11 * BLOCK];
     1662      t3 -= region[11] * a[3 + 11 * BLOCK];
     1663
     1664      t0 -= region[12] * a[0 + 12 * BLOCK];
     1665      t1 -= region[12] * a[1 + 12 * BLOCK];
     1666      t2 -= region[12] * a[2 + 12 * BLOCK];
     1667      t3 -= region[12] * a[3 + 12 * BLOCK];
     1668
     1669      t0 -= region[13] * a[0 + 13 * BLOCK];
     1670      t1 -= region[13] * a[1 + 13 * BLOCK];
     1671      t2 -= region[13] * a[2 + 13 * BLOCK];
     1672      t3 -= region[13] * a[3 + 13 * BLOCK];
     1673
     1674      t0 -= region[14] * a[0 + 14 * BLOCK];
     1675      t1 -= region[14] * a[1 + 14 * BLOCK];
     1676      t2 -= region[14] * a[2 + 14 * BLOCK];
     1677      t3 -= region[14] * a[3 + 14 * BLOCK];
     1678
     1679      t0 -= region[15] * a[0 + 15 * BLOCK];
     1680      t1 -= region[15] * a[1 + 15 * BLOCK];
     1681      t2 -= region[15] * a[2 + 15 * BLOCK];
     1682      t3 -= region[15] * a[3 + 15 * BLOCK];
     1683#endif
     1684      region2[0] = t0;
     1685      region2[1] = t1;
     1686      region2[2] = t2;
     1687      region2[3] = t3;
     1688      region2+=4;
     1689      a+=4;
     1690    }
     1691  } else {
     1692#endif
     1693    for (k = 0; k < n; ++k) {
     1694      longDouble t00 = region2[k];
     1695      for (j = 0; j < BLOCK; j ++) {
     1696        t00 -= region[j] * a[k + j * BLOCK];
     1697      }
     1698      region2[k] = t00;
     1699    }
     1700#if BLOCK==16
     1701  }
     1702#endif
     1703}
     1704// Backward part of solve 1
     1705void
     1706ClpCholeskyDense::solveB1Long(longDouble * a,int n,longDouble * region)
     1707{
     1708  int j, k;
     1709  longDouble t00;
     1710  for (j = n-1; j >=0; j --) {
     1711    t00 = region[j];
     1712    for (k = j+1; k < n; ++k) {
     1713      t00 -= region[k] * a[k + j * BLOCK];
     1714    }
     1715    //t00*=a[j + j * BLOCK];
     1716    region[j] = t00;
     1717  }
     1718}
     1719// Backward part of solve 2
     1720void
     1721ClpCholeskyDense::solveB2Long(longDouble * a,int n,longDouble * region, longDouble * region2)
     1722{
     1723  int j, k;
     1724#if BLOCK==16
     1725  if (n==BLOCK) {
     1726    for (j = 0; j < BLOCK; j +=4) {
     1727      longDouble t0 = region[0];
     1728      longDouble t1 = region[1];
     1729      longDouble t2 = region[2];
     1730      longDouble t3 = region[3];
     1731      t0 -= region2[0] * a[0 + 0*BLOCK];
     1732      t1 -= region2[0] * a[0 + 1*BLOCK];
     1733      t2 -= region2[0] * a[0 + 2*BLOCK];
     1734      t3 -= region2[0] * a[0 + 3*BLOCK];
     1735
     1736      t0 -= region2[1] * a[1 + 0*BLOCK];
     1737      t1 -= region2[1] * a[1 + 1*BLOCK];
     1738      t2 -= region2[1] * a[1 + 2*BLOCK];
     1739      t3 -= region2[1] * a[1 + 3*BLOCK];
     1740
     1741      t0 -= region2[2] * a[2 + 0*BLOCK];
     1742      t1 -= region2[2] * a[2 + 1*BLOCK];
     1743      t2 -= region2[2] * a[2 + 2*BLOCK];
     1744      t3 -= region2[2] * a[2 + 3*BLOCK];
     1745
     1746      t0 -= region2[3] * a[3 + 0*BLOCK];
     1747      t1 -= region2[3] * a[3 + 1*BLOCK];
     1748      t2 -= region2[3] * a[3 + 2*BLOCK];
     1749      t3 -= region2[3] * a[3 + 3*BLOCK];
     1750
     1751      t0 -= region2[4] * a[4 + 0*BLOCK];
     1752      t1 -= region2[4] * a[4 + 1*BLOCK];
     1753      t2 -= region2[4] * a[4 + 2*BLOCK];
     1754      t3 -= region2[4] * a[4 + 3*BLOCK];
     1755
     1756      t0 -= region2[5] * a[5 + 0*BLOCK];
     1757      t1 -= region2[5] * a[5 + 1*BLOCK];
     1758      t2 -= region2[5] * a[5 + 2*BLOCK];
     1759      t3 -= region2[5] * a[5 + 3*BLOCK];
     1760
     1761      t0 -= region2[6] * a[6 + 0*BLOCK];
     1762      t1 -= region2[6] * a[6 + 1*BLOCK];
     1763      t2 -= region2[6] * a[6 + 2*BLOCK];
     1764      t3 -= region2[6] * a[6 + 3*BLOCK];
     1765
     1766      t0 -= region2[7] * a[7 + 0*BLOCK];
     1767      t1 -= region2[7] * a[7 + 1*BLOCK];
     1768      t2 -= region2[7] * a[7 + 2*BLOCK];
     1769      t3 -= region2[7] * a[7 + 3*BLOCK];
     1770#if BLOCK>8
     1771
     1772      t0 -= region2[8] * a[8 + 0*BLOCK];
     1773      t1 -= region2[8] * a[8 + 1*BLOCK];
     1774      t2 -= region2[8] * a[8 + 2*BLOCK];
     1775      t3 -= region2[8] * a[8 + 3*BLOCK];
     1776
     1777      t0 -= region2[9] * a[9 + 0*BLOCK];
     1778      t1 -= region2[9] * a[9 + 1*BLOCK];
     1779      t2 -= region2[9] * a[9 + 2*BLOCK];
     1780      t3 -= region2[9] * a[9 + 3*BLOCK];
     1781
     1782      t0 -= region2[10] * a[10 + 0*BLOCK];
     1783      t1 -= region2[10] * a[10 + 1*BLOCK];
     1784      t2 -= region2[10] * a[10 + 2*BLOCK];
     1785      t3 -= region2[10] * a[10 + 3*BLOCK];
     1786
     1787      t0 -= region2[11] * a[11 + 0*BLOCK];
     1788      t1 -= region2[11] * a[11 + 1*BLOCK];
     1789      t2 -= region2[11] * a[11 + 2*BLOCK];
     1790      t3 -= region2[11] * a[11 + 3*BLOCK];
     1791
     1792      t0 -= region2[12] * a[12 + 0*BLOCK];
     1793      t1 -= region2[12] * a[12 + 1*BLOCK];
     1794      t2 -= region2[12] * a[12 + 2*BLOCK];
     1795      t3 -= region2[12] * a[12 + 3*BLOCK];
     1796
     1797      t0 -= region2[13] * a[13 + 0*BLOCK];
     1798      t1 -= region2[13] * a[13 + 1*BLOCK];
     1799      t2 -= region2[13] * a[13 + 2*BLOCK];
     1800      t3 -= region2[13] * a[13 + 3*BLOCK];
     1801
     1802      t0 -= region2[14] * a[14 + 0*BLOCK];
     1803      t1 -= region2[14] * a[14 + 1*BLOCK];
     1804      t2 -= region2[14] * a[14 + 2*BLOCK];
     1805      t3 -= region2[14] * a[14 + 3*BLOCK];
     1806
     1807      t0 -= region2[15] * a[15 + 0*BLOCK];
     1808      t1 -= region2[15] * a[15 + 1*BLOCK];
     1809      t2 -= region2[15] * a[15 + 2*BLOCK];
     1810      t3 -= region2[15] * a[15 + 3*BLOCK];
     1811#endif
     1812      region[0] = t0;
     1813      region[1] = t1;
     1814      region[2] = t2;
     1815      region[3] = t3;
     1816      a+=4*BLOCK;
     1817      region+=4;
     1818    }
     1819  } else {
     1820#endif
     1821    for (j = 0; j < BLOCK; j ++) {
     1822      longDouble t00 = region[j];
     1823      for (k = 0; k < n; ++k) {
     1824        t00 -= region2[k] * a[k + j * BLOCK];
     1825      }
     1826      region[j] = t00;
     1827    }
     1828#if BLOCK==16
     1829  }
     1830#endif
     1831}
     1832/* Uses factorization to solve. */
     1833void
     1834ClpCholeskyDense::solveLongWork (longWork * region)
     1835{
     1836  int numberBlocks = (numberRows_+BLOCK-1)>>BLOCKSHIFT;
     1837  // later align on boundary
     1838  longDouble * a = sparseFactor_+BLOCKSQ*numberBlocks;
     1839  int iBlock;
     1840  longDouble * aa = a;
     1841  int iColumn;
     1842  for (iBlock=0;iBlock<numberBlocks;iBlock++) {
     1843    int nChunk;
     1844    int jBlock;
     1845    int iDo = iBlock*BLOCK;
     1846    int base=iDo;
     1847    if (iDo+BLOCK>numberRows_) {
     1848      nChunk=numberRows_-iDo;
     1849    } else {
     1850      nChunk=BLOCK;
     1851    }
     1852    solveF1LongWork(aa,nChunk,region+iDo);
     1853    for (jBlock=iBlock+1;jBlock<numberBlocks;jBlock++) {
     1854      base+=BLOCK;
     1855      aa+=BLOCKSQ;
     1856      if (base+BLOCK>numberRows_) {
     1857        nChunk=numberRows_-base;
     1858      } else {
     1859        nChunk=BLOCK;
     1860      }
     1861      solveF2LongWork(aa,nChunk,region+iDo,region+base);
     1862    }
     1863    aa+=BLOCKSQ;
     1864  }
     1865  // do diagonal outside
     1866  for (iColumn=0;iColumn<numberRows_;iColumn++)
     1867    region[iColumn] *= diagonal_[iColumn];
     1868  int offset=((numberBlocks*(numberBlocks+1))>>1);
     1869  aa = a+number_entries(offset-1);
     1870  int lBase=(numberBlocks-1)*BLOCK;
     1871  for (iBlock=numberBlocks-1;iBlock>=0;iBlock--) {
     1872    int nChunk;
     1873    int jBlock;
     1874    int triBase=iBlock*BLOCK;
     1875    int iBase=lBase;
     1876    for (jBlock=iBlock+1;jBlock<numberBlocks;jBlock++) {
     1877      if (iBase+BLOCK>numberRows_) {
     1878        nChunk=numberRows_-iBase;
     1879      } else {
     1880        nChunk=BLOCK;
     1881      }
     1882      solveB2LongWork(aa,nChunk,region+triBase,region+iBase);
     1883      iBase-=BLOCK;
     1884      aa-=BLOCKSQ;
     1885    }
     1886    if (triBase+BLOCK>numberRows_) {
     1887      nChunk=numberRows_-triBase;
     1888    } else {
     1889      nChunk=BLOCK;
     1890    }
     1891    solveB1LongWork(aa,nChunk,region+triBase);
     1892    aa-=BLOCKSQ;
     1893  }
     1894}
     1895// Forward part of solve 1
     1896void
     1897ClpCholeskyDense::solveF1LongWork(longDouble * a,int n,longWork * region)
     1898{
     1899  int j, k;
     1900  longWork t00;
     1901  for (j = 0; j < n; j ++) {
     1902    t00 = region[j];
     1903    for (k = 0; k < j; ++k) {
     1904      t00 -= region[k] * a[j + k * BLOCK];
     1905    }
     1906    //t00*=a[j + j * BLOCK];
     1907    region[j] = t00;
     1908  }
     1909}
     1910// Forward part of solve 2
     1911void
     1912ClpCholeskyDense::solveF2LongWork(longDouble * a,int n,longWork * region, longWork * region2)
     1913{
     1914  int j, k;
     1915#if BLOCK==16
     1916  if (n==BLOCK) {
     1917    for (k = 0; k < BLOCK; k+=4) {
     1918      longWork t0 = region2[0];
     1919      longWork t1 = region2[1];
     1920      longWork t2 = region2[2];
     1921      longWork t3 = region2[3];
     1922      t0 -= region[0] * a[0 + 0 * BLOCK];
     1923      t1 -= region[0] * a[1 + 0 * BLOCK];
     1924      t2 -= region[0] * a[2 + 0 * BLOCK];
     1925      t3 -= region[0] * a[3 + 0 * BLOCK];
     1926
     1927      t0 -= region[1] * a[0 + 1 * BLOCK];
     1928      t1 -= region[1] * a[1 + 1 * BLOCK];
     1929      t2 -= region[1] * a[2 + 1 * BLOCK];
     1930      t3 -= region[1] * a[3 + 1 * BLOCK];
     1931
     1932      t0 -= region[2] * a[0 + 2 * BLOCK];
     1933      t1 -= region[2] * a[1 + 2 * BLOCK];
     1934      t2 -= region[2] * a[2 + 2 * BLOCK];
     1935      t3 -= region[2] * a[3 + 2 * BLOCK];
     1936
     1937      t0 -= region[3] * a[0 + 3 * BLOCK];
     1938      t1 -= region[3] * a[1 + 3 * BLOCK];
     1939      t2 -= region[3] * a[2 + 3 * BLOCK];
     1940      t3 -= region[3] * a[3 + 3 * BLOCK];
     1941
     1942      t0 -= region[4] * a[0 + 4 * BLOCK];
     1943      t1 -= region[4] * a[1 + 4 * BLOCK];
     1944      t2 -= region[4] * a[2 + 4 * BLOCK];
     1945      t3 -= region[4] * a[3 + 4 * BLOCK];
     1946
     1947      t0 -= region[5] * a[0 + 5 * BLOCK];
     1948      t1 -= region[5] * a[1 + 5 * BLOCK];
     1949      t2 -= region[5] * a[2 + 5 * BLOCK];
     1950      t3 -= region[5] * a[3 + 5 * BLOCK];
     1951
     1952      t0 -= region[6] * a[0 + 6 * BLOCK];
     1953      t1 -= region[6] * a[1 + 6 * BLOCK];
     1954      t2 -= region[6] * a[2 + 6 * BLOCK];
     1955      t3 -= region[6] * a[3 + 6 * BLOCK];
     1956
     1957      t0 -= region[7] * a[0 + 7 * BLOCK];
     1958      t1 -= region[7] * a[1 + 7 * BLOCK];
     1959      t2 -= region[7] * a[2 + 7 * BLOCK];
     1960      t3 -= region[7] * a[3 + 7 * BLOCK];
     1961#if BLOCK>8
     1962      t0 -= region[8] * a[0 + 8 * BLOCK];
     1963      t1 -= region[8] * a[1 + 8 * BLOCK];
     1964      t2 -= region[8] * a[2 + 8 * BLOCK];
     1965      t3 -= region[8] * a[3 + 8 * BLOCK];
     1966
     1967      t0 -= region[9] * a[0 + 9 * BLOCK];
     1968      t1 -= region[9] * a[1 + 9 * BLOCK];
     1969      t2 -= region[9] * a[2 + 9 * BLOCK];
     1970      t3 -= region[9] * a[3 + 9 * BLOCK];
     1971
     1972      t0 -= region[10] * a[0 + 10 * BLOCK];
     1973      t1 -= region[10] * a[1 + 10 * BLOCK];
     1974      t2 -= region[10] * a[2 + 10 * BLOCK];
     1975      t3 -= region[10] * a[3 + 10 * BLOCK];
     1976
     1977      t0 -= region[11] * a[0 + 11 * BLOCK];
     1978      t1 -= region[11] * a[1 + 11 * BLOCK];
     1979      t2 -= region[11] * a[2 + 11 * BLOCK];
     1980      t3 -= region[11] * a[3 + 11 * BLOCK];
     1981
     1982      t0 -= region[12] * a[0 + 12 * BLOCK];
     1983      t1 -= region[12] * a[1 + 12 * BLOCK];
     1984      t2 -= region[12] * a[2 + 12 * BLOCK];
     1985      t3 -= region[12] * a[3 + 12 * BLOCK];
     1986
     1987      t0 -= region[13] * a[0 + 13 * BLOCK];
     1988      t1 -= region[13] * a[1 + 13 * BLOCK];
     1989      t2 -= region[13] * a[2 + 13 * BLOCK];
     1990      t3 -= region[13] * a[3 + 13 * BLOCK];
     1991
     1992      t0 -= region[14] * a[0 + 14 * BLOCK];
     1993      t1 -= region[14] * a[1 + 14 * BLOCK];
     1994      t2 -= region[14] * a[2 + 14 * BLOCK];
     1995      t3 -= region[14] * a[3 + 14 * BLOCK];
     1996
     1997      t0 -= region[15] * a[0 + 15 * BLOCK];
     1998      t1 -= region[15] * a[1 + 15 * BLOCK];
     1999      t2 -= region[15] * a[2 + 15 * BLOCK];
     2000      t3 -= region[15] * a[3 + 15 * BLOCK];
     2001#endif
     2002      region2[0] = t0;
     2003      region2[1] = t1;
     2004      region2[2] = t2;
     2005      region2[3] = t3;
     2006      region2+=4;
     2007      a+=4;
     2008    }
     2009  } else {
     2010#endif
     2011    for (k = 0; k < n; ++k) {
     2012      longWork t00 = region2[k];
     2013      for (j = 0; j < BLOCK; j ++) {
     2014        t00 -= region[j] * a[k + j * BLOCK];
     2015      }
     2016      region2[k] = t00;
     2017    }
     2018#if BLOCK==16
     2019  }
     2020#endif
     2021}
     2022// Backward part of solve 1
     2023void
     2024ClpCholeskyDense::solveB1LongWork(longDouble * a,int n,longWork * region)
     2025{
     2026  int j, k;
     2027  longWork t00;
     2028  for (j = n-1; j >=0; j --) {
     2029    t00 = region[j];
     2030    for (k = j+1; k < n; ++k) {
     2031      t00 -= region[k] * a[k + j * BLOCK];
     2032    }
     2033    //t00*=a[j + j * BLOCK];
     2034    region[j] = t00;
     2035  }
     2036}
     2037// Backward part of solve 2
     2038void
     2039ClpCholeskyDense::solveB2LongWork(longDouble * a,int n,longWork * region, longWork * region2)
     2040{
     2041  int j, k;
     2042#if BLOCK==16
     2043  if (n==BLOCK) {
     2044    for (j = 0; j < BLOCK; j +=4) {
     2045      longWork t0 = region[0];
     2046      longWork t1 = region[1];
     2047      longWork t2 = region[2];
     2048      longWork t3 = region[3];
     2049      t0 -= region2[0] * a[0 + 0*BLOCK];
     2050      t1 -= region2[0] * a[0 + 1*BLOCK];
     2051      t2 -= region2[0] * a[0 + 2*BLOCK];
     2052      t3 -= region2[0] * a[0 + 3*BLOCK];
     2053
     2054      t0 -= region2[1] * a[1 + 0*BLOCK];
     2055      t1 -= region2[1] * a[1 + 1*BLOCK];
     2056      t2 -= region2[1] * a[1 + 2*BLOCK];
     2057      t3 -= region2[1] * a[1 + 3*BLOCK];
     2058
     2059      t0 -= region2[2] * a[2 + 0*BLOCK];
     2060      t1 -= region2[2] * a[2 + 1*BLOCK];
     2061      t2 -= region2[2] * a[2 + 2*BLOCK];
     2062      t3 -= region2[2] * a[2 + 3*BLOCK];
     2063
     2064      t0 -= region2[3] * a[3 + 0*BLOCK];
     2065      t1 -= region2[3] * a[3 + 1*BLOCK];
     2066      t2 -= region2[3] * a[3 + 2*BLOCK];
     2067      t3 -= region2[3] * a[3 + 3*BLOCK];
     2068
     2069      t0 -= region2[4] * a[4 + 0*BLOCK];
     2070      t1 -= region2[4] * a[4 + 1*BLOCK];
     2071      t2 -= region2[4] * a[4 + 2*BLOCK];
     2072      t3 -= region2[4] * a[4 + 3*BLOCK];
     2073
     2074      t0 -= region2[5] * a[5 + 0*BLOCK];
     2075      t1 -= region2[5] * a[5 + 1*BLOCK];
     2076      t2 -= region2[5] * a[5 + 2*BLOCK];
     2077      t3 -= region2[5] * a[5 + 3*BLOCK];
     2078
     2079      t0 -= region2[6] * a[6 + 0*BLOCK];
     2080      t1 -= region2[6] * a[6 + 1*BLOCK];
     2081      t2 -= region2[6] * a[6 + 2*BLOCK];
     2082      t3 -= region2[6] * a[6 + 3*BLOCK];
     2083
     2084      t0 -= region2[7] * a[7 + 0*BLOCK];
     2085      t1 -= region2[7] * a[7 + 1*BLOCK];
     2086      t2 -= region2[7] * a[7 + 2*BLOCK];
     2087      t3 -= region2[7] * a[7 + 3*BLOCK];
     2088#if BLOCK>8
     2089
     2090      t0 -= region2[8] * a[8 + 0*BLOCK];
     2091      t1 -= region2[8] * a[8 + 1*BLOCK];
     2092      t2 -= region2[8] * a[8 + 2*BLOCK];
     2093      t3 -= region2[8] * a[8 + 3*BLOCK];
     2094
     2095      t0 -= region2[9] * a[9 + 0*BLOCK];
     2096      t1 -= region2[9] * a[9 + 1*BLOCK];
     2097      t2 -= region2[9] * a[9 + 2*BLOCK];
     2098      t3 -= region2[9] * a[9 + 3*BLOCK];
     2099
     2100      t0 -= region2[10] * a[10 + 0*BLOCK];
     2101      t1 -= region2[10] * a[10 + 1*BLOCK];
     2102      t2 -= region2[10] * a[10 + 2*BLOCK];
     2103      t3 -= region2[10] * a[10 + 3*BLOCK];
     2104
     2105      t0 -= region2[11] * a[11 + 0*BLOCK];
     2106      t1 -= region2[11] * a[11 + 1*BLOCK];
     2107      t2 -= region2[11] * a[11 + 2*BLOCK];
     2108      t3 -= region2[11] * a[11 + 3*BLOCK];
     2109
     2110      t0 -= region2[12] * a[12 + 0*BLOCK];
     2111      t1 -= region2[12] * a[12 + 1*BLOCK];
     2112      t2 -= region2[12] * a[12 + 2*BLOCK];
     2113      t3 -= region2[12] * a[12 + 3*BLOCK];
     2114
     2115      t0 -= region2[13] * a[13 + 0*BLOCK];
     2116      t1 -= region2[13] * a[13 + 1*BLOCK];
     2117      t2 -= region2[13] * a[13 + 2*BLOCK];
     2118      t3 -= region2[13] * a[13 + 3*BLOCK];
     2119
     2120      t0 -= region2[14] * a[14 + 0*BLOCK];
     2121      t1 -= region2[14] * a[14 + 1*BLOCK];
     2122      t2 -= region2[14] * a[14 + 2*BLOCK];
     2123      t3 -= region2[14] * a[14 + 3*BLOCK];
     2124
     2125      t0 -= region2[15] * a[15 + 0*BLOCK];
     2126      t1 -= region2[15] * a[15 + 1*BLOCK];
     2127      t2 -= region2[15] * a[15 + 2*BLOCK];
     2128      t3 -= region2[15] * a[15 + 3*BLOCK];
     2129#endif
     2130      region[0] = t0;
     2131      region[1] = t1;
     2132      region[2] = t2;
     2133      region[3] = t3;
     2134      a+=4*BLOCK;
     2135      region+=4;
     2136    }
     2137  } else {
     2138#endif
     2139    for (j = 0; j < BLOCK; j ++) {
     2140      longWork t00 = region[j];
     2141      for (k = 0; k < n; ++k) {
     2142        t00 -= region2[k] * a[k + j * BLOCK];
     2143      }
     2144      region[j] = t00;
     2145    }
     2146#if BLOCK==16
     2147  }
     2148#endif
     2149}
  • trunk/ClpCholeskyTaucs.cpp

    r414 r423  
    1 #ifdef REAL_BARRIER
     1#ifdef TAUCS_BARRIER
    22// Copyright (C) 2004, International Business Machines
    33// Corporation and others.  All Rights Reserved.
     
    230230  taucs_ccs_free(permuted);
    231231  return factorization_ ? 0 :1;
     232}
     233/* Does Symbolic factorization given permutation.
     234   This is called immediately after order.  If user provides this then
     235   user must provide factorize and solve.  Otherwise the default factorization is used
     236   returns non-zero if not enough memory */
     237int
     238ClpCholeskyTaucs::symbolic()
     239{
     240  return 0;
    232241}
    233242/* Factorize - filling in rowsDropped and returning number dropped */
  • trunk/ClpCholeskyWssmp.cpp

    r414 r423  
    1 #ifdef REAL_BARRIER
     1#ifdef WSSMP_BARRIER
    22// Copyright (C) 2003, International Business Machines
    33// Corporation and others.  All Rights Reserved.
     
    6363// At present I can't get wssmp to work as my libraries seem to be out of sync
    6464// so I have linked in ekkwssmp which is an older version
    65 //#define WSMP
    66 #ifdef WSMP
     65#if WSSMP_BARRIER==2
    6766  extern "C" void wssmp(int * n,
    6867                        int * columnStart , int * rowIndex , double * element,
     
    277276  int i0=0;
    278277  int i1=1;
    279 #ifdef WSMP   
     278#if WSSMP_BARRIER==2
    280279  wsetmaxthrds(&i1);
    281280#endif
  • trunk/ClpCholeskyWssmpKKT.cpp

    r414 r423  
    1 #ifdef REAL_BARRIER
     1#ifdef WSSMP_BARRIER
    22// Copyright (C) 2004, International Business Machines
    33// Corporation and others.  All Rights Reserved.
     
    6363// At present I can't get wssmp to work as my libraries seem to be out of sync
    6464// so I have linked in ekkwssmp which is an older version
    65 //#define WSMP
    66 #ifdef WSMP
     65#if WSSMP_BARRIER==2
    6766  extern "C" void wssmp(int * n,
    6867                        int * columnStart , int * rowIndex , double * element,
     
    212211  int i0=0;
    213212  int i1=1;
    214 #ifdef WSMP   
     213#if WSSMP_BARRIER==2
    215214  wsetmaxthrds(&i1);
    216215#endif
  • trunk/ClpSimplex.cpp

    r414 r423  
    35663566}
    35673567#include "ClpPredictorCorrector.hpp"
    3568 #ifdef REAL_BARRIER
     3568#include "ClpCholeskyBase.hpp"
     3569// Preference is WSSMP, TAUCS, UFL (just ordering) then base
     3570#if WSSMP_BARRIER
    35693571#include "ClpCholeskyWssmp.hpp"
    35703572#include "ClpCholeskyWssmpKKT.hpp"
    3571 //#include "ClpCholeskyTaucs.hpp"
     3573#elif UFL_BARRIER
     3574#include "ClpCholeskyUfl.hpp"
     3575#elif TAUCS_BARRIER
     3576#include "ClpCholeskyTaucs.hpp"
    35723577#endif
    35733578#include "ClpPresolve.hpp"
     
    35813586  ClpInterior barrier;
    35823587  barrier.borrowModel(*model2);
    3583 #ifdef REAL_BARRIER
    3584   // uncomment this if you have Anshul Gupta's wsmp package
     3588// Preference is WSSMP, UFL, TAUCS then base
     3589#ifdef WSSMP_BARRIER
    35853590  ClpCholeskyWssmp * cholesky = new ClpCholeskyWssmp(CoinMax(100,model2->numberRows()/10));
    35863591  //ClpCholeskyWssmp * cholesky = new ClpCholeskyWssmp();
    35873592  //ClpCholeskyWssmpKKT * cholesky = new ClpCholeskyWssmpKKT(CoinMax(100,model2->numberRows()/10));
    3588   // uncomment this if you have Sivan Toledo's Taucs package
    3589   //ClpCholeskyTaucs * cholesky = new ClpCholeskyTaucs();
     3593#elif UFL_BARRIER
     3594  ClpCholeskyUfl * cholesky = new ClpCholeskyUfl();
     3595#elif TAUCS_BARRIER
     3596  ClpCholeskyTaucs * cholesky = new ClpCholeskyTaucs();
     3597#else
     3598  ClpCholeskyBase * cholesky = new ClpCholeskyBase();
     3599#endif
    35903600  barrier.setCholesky(cholesky);
    3591 #endif
    35923601  int numberRows = model2->numberRows();
    35933602  int numberColumns = model2->numberColumns();
  • trunk/ClpSolve.cpp

    r414 r423  
    1818#include "ClpCholeskyDense.hpp"
    1919#include "ClpCholeskyBase.hpp"
    20 #ifdef REAL_BARRIER
    21 #include "ClpCholeskyWssmp.hpp"
    22 #include "ClpCholeskyWssmpKKT.hpp"
    23 //#include "ClpCholeskyTaucs.hpp"
    24 #include "ClpCholeskyUfl.hpp"
    25 #else
    26 static int numberBarrier=0;
    27 #endif
    2820#include "ClpSolve.hpp"
    2921#include "ClpPackedMatrix.hpp"
     
    3527#include "ClpPresolve.hpp"
    3628#include "Idiot.hpp"
     29#ifdef WSSMP_BARRIER
     30#include "ClpCholeskyWssmp.hpp"
     31#include "ClpCholeskyWssmpKKT.hpp"
     32#define FAST_BARRIER
     33#endif
     34#ifdef UFL_BARRIER
     35#include "ClpCholeskyUfl.hpp"
     36#endif
     37#ifdef TAUCS_BARRIER
     38#include "ClpCholeskyTaucs.hpp"
     39#define FAST_BARRIER
     40#endif
     41#ifndef FAST_BARRIER
     42static int numberBarrier=0;
     43#endif
    3744
    3845//#############################################################################
     
    949956      quadraticObj = (static_cast< ClpQuadraticObjective*>(objective_));
    950957#endif
    951 #ifdef REAL_BARRIER
     958#ifndef FAST_BARRIER
     959    if (!numberBarrier)
     960      std::cout<<"Warning - the default ordering is just on row counts! "
     961               <<"The factorization is being improved"<<std::endl;
     962    numberBarrier++;
     963#endif
    952964    if (quadraticObj) {
    953965      doKKT=true;
     
    955967    switch (barrierOptions) {
    956968    case 0:
     969    default:
    957970      if (!doKKT) {
    958         ClpCholeskyBase * cholesky = new ClpCholeskyBase(100);
     971        ClpCholeskyBase * cholesky = new ClpCholeskyBase();
    959972        barrier.setCholesky(cholesky);
    960973      } else {
     
    974987      }
    975988      break;
     989#ifdef WSSMP_BARRIER
    976990    case 2:
    977991      {
     
    9901004      }
    9911005      break;
     1006#endif
     1007#ifdef UFL_BARRIER
    9921008    case 4:
    9931009      if (!doKKT) {
    994         ClpCholeskyUfl * cholesky = new ClpCholeskyUfl(100);
     1010        ClpCholeskyUfl * cholesky = new ClpCholeskyUfl();
    9951011        barrier.setCholesky(cholesky);
    9961012      } else {
     
    10001016      }
    10011017      break;
    1002     default:
    1003       abort();
    1004     }
    1005     // uncomment this if you have Sivan Toledo's Taucs package
    1006     //ClpCholeskyTaucs * cholesky = new ClpCholeskyTaucs();
    1007     //barrier.setCholesky(cholesky);
    1008 #else
    1009     if (quadraticObj) {
    1010       doKKT=true;
    1011     }
    1012     switch (barrierOptions) {
    1013     case 1:
    1014       if (!doKKT) {
    1015         ClpCholeskyDense * cholesky = new ClpCholeskyDense();
     1018#endif
     1019#ifdef TAUCS_BARRIER
     1020    case 5:
     1021      {
     1022        ClpCholeskyTaucs * cholesky = new ClpCholeskyTaucs();
    10161023        barrier.setCholesky(cholesky);
    1017       } else {
    1018         ClpCholeskyDense * cholesky = new ClpCholeskyDense();
    1019         cholesky->setKKT(true);
    1020         barrier.setCholesky(cholesky);
    1021       }
    1022       break;
    1023     default:
    1024       if (!numberBarrier)
    1025         std::cout<<"Warning - the default ordering is just on row counts! "
    1026                  <<"The factorization is being improved"<<std::endl;
    1027       numberBarrier++;
    1028       if (!doKKT) {
    1029         ClpCholeskyBase * cholesky = new ClpCholeskyBase(100);
    1030         barrier.setCholesky(cholesky);
    1031       } else {
    1032         ClpCholeskyBase * cholesky = new ClpCholeskyBase();
    1033         cholesky->setKKT(true);
    1034         barrier.setCholesky(cholesky);
    1035       }
    1036       break;
    1037     }
    1038 #endif
     1024        assert (!doKKT);
     1025      }
     1026      break;
     1027#endif
     1028    }
    10391029    int numberRows = model2->numberRows();
    10401030    int numberColumns = model2->numberColumns();
  • trunk/Samples/decompose.cpp

    r343 r423  
    1313  if (argc<2) {
    1414    //status=model.readMps("/home/forrest/data/ken-18.mps.gz",true);
    15     status=model.readMps("small.mps",true);
     15    status=model.readMps("../../Mps/Netlib/czprob.mps",true);
    1616  } else {
    1717    status=model.readMps(argv[1],true);
     
    5050    // degen3
    5151    for (iRow=0;iRow<561;iRow++)
     52      rowBlock[iRow]=-1;
     53  } else if (numberRows==929) {
     54    // czprob
     55    for (iRow=0;iRow<39;iRow++)
    5256      rowBlock[iRow]=-1;
    5357  }
     
    184188  numberColumn2=0;
    185189  for (iRow=0;iRow<numberRows;iRow++)
    186     if (rowBlock[iRow]==-1)
     190    if (rowBlock[iRow]<0)
    187191      whichRow[numberRow2++]=iRow;
    188192  for (iColumn=0;iColumn<numberColumns;iColumn++)
  • trunk/Test/ClpMain.cpp

    r414 r423  
    4040#ifdef DMALLOC
    4141#include "dmalloc.h"
     42#endif
     43#ifdef WSSMP_BARRIER
     44#define FOREIGN_BARRIER
     45#endif
     46#ifdef UFL_BARRIER
     47#define FOREIGN_BARRIER
     48#endif
     49#ifdef TAUCS_BARRIER
     50#define FOREIGN_BARRIER
    4251#endif
    4352
     
    942951    parameters[numberParameters++]=
    943952      ClpItem("chol!esky","Which cholesky algorithm",
    944               "base",CHOLESKY,false);
     953              "native",CHOLESKY,false);
    945954    parameters[numberParameters-1].append("dense");
    946 #ifdef REAL_BARRIER
     955#ifdef FOREIGN_BARRIER
     956#ifdef WSSMP_BARRIER
    947957    parameters[numberParameters-1].append("fudge!Long");
    948958    parameters[numberParameters-1].append("wssmp");
     959#else
     960    parameters[numberParameters-1].append("fudge!Long_dummy");
     961    parameters[numberParameters-1].append("wssmp_dummy");
     962#endif
     963#ifdef UFL_BARRIER
    949964    parameters[numberParameters-1].append("Uni!versityOfFlorida");
     965#else
     966    parameters[numberParameters-1].append("Uni!versityOfFlorida_dummy");
     967#endif
     968#ifdef TAUCS_BARRIER
     969    parameters[numberParameters-1].append("Taucs");
     970#else
     971    parameters[numberParameters-1].append("Taucs_dummy");
     972#endif
    950973#endif
    951974    parameters[numberParameters++]=
  • trunk/include/ClpCholeskyBase.hpp

    r414 r423  
    66#include "CoinPragma.hpp"
    77#include "CoinFinite.hpp"
    8 
     8#define CLP_LONG_CHOLESKY 0
     9#if CLP_LONG_CHOLESKY>1
     10typedef long double longDouble;
     11typedef long double longWork;
     12#elif CLP_LONG_CHOLESKY==1
     13typedef double longDouble;
     14typedef long double longWork;
     15#else
     16typedef double longDouble;
     17typedef double longWork;
     18#endif
    919class ClpInterior;
    1020
     
    1222    Will do better factorization.  very crude ordering
    1323
    14     Derived classes will be using sophisticated methods apart from ClpCholeskyDense
     24    Derived classes may be using more sophisticated methods
    1525*/
    1626
     
    6171  inline double choleskyCondition() const
    6272  {return choleskyCondition_;};
     73  /// goDense i.e. use dense factoriaztion if > this (default 0.7).
     74  inline double goDense() const
     75  {return goDense_;};
     76  /// goDense i.e. use dense factoriaztion if > this (default 0.7).
     77  inline void setGoDense(double value)
     78  {goDense_=value;};
    6379  /// rank.  Returns rank
    6480  inline int rank() const
     
    7187  { return sizeFactor_;};
    7288  /// Return sparseFactor
    73   inline double * sparseFactor() const
     89  inline longDouble * sparseFactor() const
    7490  { return sparseFactor_;};
    7591  /// Return diagonal
    76   inline double * diagonal() const
     92  inline longDouble * diagonal() const
    7793  { return diagonal_;};
    7894  /// Return workDouble
    79   inline double * workDouble() const
     95  inline longDouble * workDouble() const
    8096  { return workDouble_;};
    8197  /// If KKT on
     
    85101  inline void setKKT(bool yesNo)
    86102  { doKKT_ = yesNo;};
     103  /// Set integer parameter
     104  inline void setIntegerParameter(int i,int value)
     105  { integerParameters_[i]=value;};
     106  /// get integer parameter
     107  inline int getIntegerParameter(int i)
     108  { return integerParameters_[i];};
     109  /// Set double parameter
     110  inline void setDoubleParameter(int i,double value)
     111  { doubleParameters_[i]=value;};
     112  /// get double parameter
     113  inline double getDoubleParameter(int i)
     114  { return doubleParameters_[i];};
    87115   //@}
    88116 
     
    135163  */
    136164  void solve(double * region, int type);
     165  void solveLong(longDouble * region, int type);
    137166  /// Forms ADAT - returns nonzero if not enough memory
    138167  int preOrder(bool lowerTriangular, bool includeDiagonal, bool doKKT);
     168  /// Updates dense part (broken out for profiling)
     169  void updateDense(longDouble * d, longDouble * work, int * first);
    139170  //@}
    140171   
     
    147178  /// Doing full KKT (only used if default symbolic and factorization)
    148179  bool doKKT_;
    149   /// zeroTolerance.
    150   double zeroTolerance_;
     180  /// Go dense at this fraction
     181  double goDense_;
    151182  /// choleskyCondition.
    152183  double choleskyCondition_;
     
    168199  int numberRowsDropped_;
    169200  /// sparseFactor.
    170   double * sparseFactor_;
     201  longDouble * sparseFactor_;
    171202  /// choleskyStart - element starts
    172203  CoinBigIndex * choleskyStart_;
     
    176207  CoinBigIndex * indexStart_;
    177208  /// Diagonal
    178   double * diagonal_;
     209  longDouble * diagonal_;
    179210  /// double work array
    180   double * workDouble_;
     211  longDouble * workDouble_;
    181212  /// link array
    182213  int * link_;
     
    200231  char * whichDense_;
    201232  /// Dense columns (updated)
    202   double * denseColumn_;
     233  longDouble * denseColumn_;
    203234  /// Dense cholesky
    204235  ClpCholeskyDense * dense_;
    205   /// Dense threshold
     236  /// Dense threshold (for taking out of Cholesky)
    206237  int denseThreshold_;
    207238  //@}
  • trunk/include/ClpCholeskyDense.hpp

    r414 r423  
    55
    66#include "ClpCholeskyBase.hpp"
    7 
    8 typedef double longDouble;
    97
    108/** Dense class for Clp Cholesky factorization
     
    4543  void factorizePart3(int * rowsDropped) ;
    4644  /// Non leaf recursive factor
    47   void factor(double * a, int n, int numberBlocks,
    48               double * diagonal, double * work, int * rowsDropped);
     45  void factor(longDouble * a, int n, int numberBlocks,
     46              longDouble * diagonal, longDouble * work, int * rowsDropped);
    4947  /// Non leaf recursive triangle rectangle update
    50   void triRec(double * aTri, int nThis, double * aUnder, double * diagonal, double * work,
     48  void triRec(longDouble * aTri, int nThis, longDouble * aUnder, longDouble * diagonal, longDouble * work,
    5149              int nLeft, int iBlock, int jBlock,
    5250              int numberBlocks);
    5351  /// Non leaf recursive rectangle triangle update
    54   void recTri(double * aUnder, int nTri, int nDo,
    55               int iBlock, int jBlock,double * aTri,
    56               double * diagonal, double * work,
     52  void recTri(longDouble * aUnder, int nTri, int nDo,
     53              int iBlock, int jBlock,longDouble * aTri,
     54              longDouble * diagonal, longDouble * work,
    5755              int numberBlocks);
    5856  /** Non leaf recursive rectangle rectangle update,
     
    6058      nUnderK is number of rows in kBlock
    6159  */
    62   void recRec(double * above, int nUnder, int nUnderK,
    63               int nDo, double * aUnder, double *aOther,double * diagonal, double * work,
     60  void recRec(longDouble * above, int nUnder, int nUnderK,
     61              int nDo, longDouble * aUnder, longDouble *aOther,longDouble * diagonal, longDouble * work,
    6462              int kBlock,int iBlock, int jBlock,
    6563              int numberBlocks);
    6664  /// Leaf recursive factor
    67   void factorLeaf(double * a, int n,
    68               double * diagonal, double * work, int * rowsDropped);
     65  void factorLeaf(longDouble * a, int n,
     66              longDouble * diagonal, longDouble * work, int * rowsDropped);
    6967  /// Leaf recursive triangle rectangle update
    70   void triRecLeaf(double * aTri, double * aUnder,
    71                   double * diagonal, double * work,
     68  void triRecLeaf(longDouble * aTri, longDouble * aUnder,
     69                  longDouble * diagonal, longDouble * work,
    7270                  int nUnder);
    7371  /// Leaf recursive rectangle triangle update
    74   void recTriLeaf(double * aUnder, double * aTri,
    75                   double * diagonal, double * work, int nUnder);
     72  void recTriLeaf(longDouble * aUnder, longDouble * aTri,
     73                  longDouble * diagonal, longDouble * work, int nUnder);
    7674  /** Leaf recursive rectangle rectangle update,
    7775      nUnder is number of rows in iBlock,
    7876      nUnderK is number of rows in kBlock
    7977  */
    80   void recRecLeaf(double * above,
    81                   double * aUnder, double *aOther, double * diagonal, double * work,
     78  void recRecLeaf(longDouble * above,
     79                  longDouble * aUnder, longDouble *aOther, longDouble * diagonal, longDouble * work,
    8280                  int nUnder);
    8381  /// Forward part of solve
    84   void solveF1(double * a,int n,double * region);
    85   void solveF2(double * a,int n,double * region,double * region2);
     82  void solveF1(longDouble * a,int n,double * region);
     83  void solveF2(longDouble * a,int n,double * region,double * region2);
    8684  /// Backward part of solve
    87   void solveB1(double * a,int n,double * region);
    88   void solveB2(double * a,int n,double * region,double * region2);
    89   int bNumber(const double * array,int &, int&);
     85  void solveB1(longDouble * a,int n,double * region);
     86  void solveB2(longDouble * a,int n,double * region,double * region2);
     87  /** Uses factorization to solve. */
     88  void solveLong (longDouble * region) ;
     89  /// Forward part of solve
     90  void solveF1Long(longDouble * a,int n,longDouble * region);
     91  void solveF2Long(longDouble * a,int n,longDouble * region,longDouble * region2);
     92  /// Backward part of solve
     93  void solveB1Long(longDouble * a,int n,longDouble * region);
     94  void solveB2Long(longDouble * a,int n,longDouble * region,longDouble * region2);
     95  /** Uses factorization to solve. */
     96  void solveLongWork (longWork * region) ;
     97  /// Forward part of solve
     98  void solveF1LongWork(longDouble * a,int n,longWork * region);
     99  void solveF2LongWork(longDouble * a,int n,longWork * region,longWork * region2);
     100  /// Backward part of solve
     101  void solveB1LongWork(longDouble * a,int n,longWork * region);
     102  void solveB2LongWork(longDouble * a,int n,longWork * region,longWork * region2);
     103  int bNumber(const longDouble * array,int &, int&);
    90104  /// A
    91   inline double * aMatrix() const
     105  inline longDouble * aMatrix() const
    92106  { return sparseFactor_;}
    93107  /// Diagonal
    94   inline double * diagonal() const
     108  inline longDouble * diagonal() const
    95109  { return diagonal_;}
    96110  //@}
  • trunk/include/ClpCholeskyTaucs.hpp

    r414 r423  
    6262   Returns non-zero if not enough memory */
    6363  virtual int order(ClpInterior * model) ;
     64  /// Dummy
     65  virtual int symbolic();
    6466  /** Factorize - filling in rowsDropped and returning number dropped.
    6567      If return code negative then out of memory */
     
    9092  /// Taucs matrix (== sparseFactor etc)
    9193  taucs_ccs_matrix * matrix_;
    92   /// Taucs factot
     94  /// Taucs factor
    9395  void * factorization_;
    9496  /// sparseFactor.
Note: See TracChangeset for help on using the changeset viewer.